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SUMMARY 


Helicopters operating in high threat areas have to fly close to the earth surface to 
minimize the risk of being detected by the adversaries. This report presents techniques for 
low altitude helicopter trajectory planning. These methods are based on optimal control 
theory and appear to be implementable onboard in realtime. Second order necessary 
conditions are obtained to provide a criterion for finding the optimal trajectory when more 
than one extremal passes through a given point. A second trajectory planning method 
incorporating a quadratic performance index is also discussed. In a later part of the thesis, 
trajectory planning problem is formulated as a differential game. The objective here is to 
synthesize optimal trajectories in the presence of an actively maneuvering adversary. 
Numerical methods for obtaining solutions to these problems are outlined. As an 
alternative to numerical method, feedback linearizing transformations are combined with the 
linear quadratic game results to synthesize explicit nonlinear feedback strategies for 
helicopter pursuit-evasion. Some of the trajectories generated from this research are 
evaluated on a six-degree-of-ffeedom helicopter simulation incorporating an advanced 
autopilot. The optimal trajectory planning methods presented here are also useful for 
autonomous land vehicle guidance. 



CHAPTER I 


INTRODUCTION 
LI Introduction 

Recent years have seen an increased interest in helicopter operations near the ground as 
evident in the literature [1-3]. In a high threat environment, helicopters have to fly close to 
the earth surface to minimize the risk of being detected by the enemy [4-5]. The objective 
here is to use terrain and surrounding objects to mask the helicopter during the mission . 

Due to data processing limitations, a hierarchical system architecture is essential for 
nap-of-the-earth flight guidance. This concept provides a natural way of decomposing a 
complex control process into simpler and more manageable components. Thus, the 
guidance functions are divided into three levels, namely, far-field, mid-field, and near-field 
[ 6 ]. 

The far-field planning task involves off-line mission planning to generate mission 
way-points and goals. Mission requirements, global threat information and vehicle 
resources on-board are taken into account. The mid-field planning function generates the 
flight route using the way-points data given by the far-field planner. High resolution digital 
map, threat information, and vehicle limitations are included in performing real-time 
guidance computations. The near-field guidance function provides a least expected 
deviation path from the mid-field nominal path due to the vehicle dynamics limitations and 
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obstacles detected by on-board sensors. The focus of this report will be on the mid-field 
route planning problem. 

Most route planning methods given in the literature [7-11] appear to use the terrain 
altitude and lateral deviations from a nominal trajectory as the performance index. These 
methods are based on the heuristic search techniques including variants of dynamic 
programming such as the A*-algorithm. All these approaches employ the discretization of 
the terrain spatial coordinates before carrying out a systematic search for optimal trajectory. 
On a rough terrain, these approaches require an enormous amount of computation and 
storage to generate sufficiently smooth trajectories [12]. 

An alternative formulation for the trajectory planning is based on Pontryagin's 
maximum principle and was first outlined in Reference 13. State equations in this 
formulation include the terrain constraint, incorporated via a coordinate transformation. 
The performance index is a linear combination of flight time and terrain altitude. The 
resulting nonlinear two-point boundary value problem is then converted to a one- 
dimensional search process by incorporating a constant of motion and employing an 
adjoint-control transformation. The solution is implementable in near real time and is 
capable of detecting situations where more than one extremal passes through a given point. 
The second-order necessary condition for this problem is studied in detail. This trajectory 
planning method automatically accomplishes known-threat avoidance and is similar to the 
classical Zermelo's navigation problem [14]. In this method, the computational algorithm 
requires the second partial derivatives of the terrain profile to generate extremals. As a 
result, the terrain profile needs to be represented by quadratic or cubic splines lattices. This 
feature can sometimes make the extremals sensitive to the error in the terrain data. 

In an alternative formulation [15], the need for second partial derivatives is eliminated 
by avoiding the coordinate transformation approach. The performance index in this 
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problem consists of a quadratic form in the terrain altitude, lateral deviation from the 
nominal trajectory, and heading angle. By changing the independent variable from time to 
down-range, the order of the problem is reduced. As in the first method, the extremals are 
obtained using the optimal control theory and necessary conditions are tested along the 
extremals. An approximate second variation test is developed for this problem using the 
WKB method [16]. A special case that can result in singular arcs in this trajectory planning 
problem is also outlined. 

So far, the route planning problem for single vehicle has been discussed. As a natural 
extension, the guidance for two or more vehicles that cooperate or compete against each 
other is considered next. This results in a differential game formulation for the trajectory 
planning problem. Since the publication of a book by Isaacs [17] on differential games in 
1965, a body of research is available on differential games with kinematic models in a 
plane. With such simple modeling, it is possible to obtain elegant results. The well-known 
homicidal chauffeur problem is an example. On the other hand, the helicopter guidance 
problem requires the use of a model in which the coefficients vary as a function of the 
vehicle position on the terrain. The method proposed in this report uses the terrain profile 
data to formulate a differential game between two helicopters. 

In conjunction with the recent theory of nonlinear transformations, Menon [18] 
showed that a class of differential games with nonlinear dynamics can be transformed into 
the well known linear quadratic pursuit-evasion game form. Compared with the previous 
derivations of pursuit-evasion guidance laws which completely ignore the dynamic 
nonlinearities in the vehicle models, the nonlinear transformation approach continuously 
compensates for the vehicle nonlinearities. In the present work, this formalism is used to 
study a helicopter pursuit-evasion game at nap-of-the-earth flight altitudes. 
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Finally, in order to verify whether the trajectories generated using various planning 
schemes discussed in the foregoing satisfy the helicopter physical constraints, these need to 
be evaluated on a detailed helicopter simulation. An advanced autopilot developed by 
Heiges [19] together with a six-degree-of-freedom helicopter simulation is used in this 
investigation. The helicopter simulation was originally developed at NASA Ames Research 
Center for the study of Air- to- Air combat [20-21]. 


1.2 Contributions of the Report 

In contrast with the existing literature, this report develops techniques for trajectory 
planning based on the Calculus of Variation. Numerical algorithms are given for the 
determination of optimal trajectories with various performance indices. Additionally, tests 
are developed for verifying the optimality of the emerging trajectories. 

Methods developed in the present research will aid in constructing an integrated 
methodology for low altitude flight guidance of helicopters. The trajectory planning 
solution is also useful for autonomous surface/underwater vehicle guidance, terrain 
following guidance for cruise missiles and aircraft [22-24] and optimal trajectory planning 
for robots. 


1.3 Organization of the Report 
This report is organized as follows: 

Chapter II gives a brief description of previous research on helicopter low-altitude 
flight trajectory planning and air-to-air combat. It was the work in this area that motivated 
the present research topic. This chapter also covers a few well-known results in optimal 
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control theory and differential game theory. Background on the nonlinear transformation 
techniques to control nonlinear systems is also presented. 

Two optimal trajectory planning schemes useful for the terrain-following/terrain- 
avoidance guidance of helicopter are presented in Chapter HI. This chapter illustrates how 
the nonlinear two-point boundary value problem can be solved using a one-dimensional 
searching method. To ensure that the extremals obtained by this approach are optimal, 
second-order necessary conditions are also developed in this chapter. 

In Chapter IV, research on the helicopter pursuit-evasion is discussed. A backward 
integration method and a nonlinear transformation method are given in this chapter. 

Chapter V discusses the implementation and test of the generated trajectories in a 
realistic six degrees of freedom helicopter simulation. The helicopter physical variables 
along the trajectories obtained from Chapter HI are examined here. 

Chapter VI evaluates the results obtained from present research. Suggestions for 
future work are also outlined. 

Finally, the appendices contain some of the analysis used in the main body of the 
report. In Appendix A, the transformation from local tangent plane to inertial coordinates is 
derived. This transformation is employed in developing the first trajectory planning 
scheme (ORP #1). Various numerical conjugate point tests and their relationships are 
discussed in Appendix B. These tests are used to verify the optimality of the synthesized 
trajectories. In Appendix C, separability of the Hamiltonian and its consequence on 
Differential Game solutions are discussed. A necessary condition for global minimum is 
given in Appendix D. 

All numerical results are obtained with VAX-1 1/750™. Contour maps are drawn by 
DISSPLA™ graphics routine. Algebraic equations are derived by symbolic program 
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MACSYMA™. Unless otherwise mentioned, British Units, i.e., pound (lb) - foot (ft) - 
second (sec), are the basic units used in this report. 
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CHAPTER II 
BACKGROUND 


2.1 Introduction 

In this chapter, previous research on the helicopter trajectory planning problem and air- 
to-air combat are reviewed. An overview of Dynamic Programming used in several of 
these research is given in Section 2.4. This section also provides an outline on optimal 
control theory. Section 2.5 provides a review of several notions involved in differential 
games. Finally, some recent results in nonlinear transformations for feedback control are 
reviewed in Section 2.6. 


2.2 Previous Research o n Helicopter Trajectory Planning 

Historically, terrain information has been used for low altitude flight guidance of deep 
penetration attack aircraft and cruise missiles. Since these vehicles fly a considerable time 
over the opponent's territory, they are vulnerable to detection by the enemy. The objective 
of low altitude flight guidance using terrain map is to minimize the influence of air defense 
threats on the mission profile [25-26]. Trajectory generated by such a guidance scheme is 
composed of a terrain-following path in the vertical plane. In the nap-of-the-earth guidance 
of helicopters, on the other hand, both vertical and lateral maneuvers are employed. 
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Reference 7 discusses the computation of vertical and lateral helicopter trajectories 
using a combination of discrete dynamic programming [27] and tree searching [28]. They 
considered the performance indices for the lateral and vertical planes as follows: 

Jl = X [wD? + (hci + Hi) 2 ] (21) 


Jv = X<h c i + Hj ) 2 
i 


( 2 . 2 ) 


where, w is the terrain-following/terrain-avoidance ratio, Di the lateral deviation from 
reference path, Hj terrain altitude at location index i, and h^ helicopter clearance altitude. 
Note that the two performance indices do not include control terms. 

Reference 1 1 developed an algorithm to generate a low altitude threat penetration 
trajectory which minimizes the performance index: 

J = (Di + CJAti ^2 3) 


Here, D, is the value of the danger array at the ith cell, Atj the transition time, and Q the 
cost of time. Danger arrays Dj depends on the vehicle position (x,y,h), and the heading 
angle X- Q is a coefficient including flight time and fuel. Decoupled vertical and lateral 
threat penetration trajectories were obtained by dynamic programming and tree search. 

Reference 29 describes a three-dimensional dynamic programming approach to 
maximize the overall probability of survival P s along any path defined by: 

P s = ri P s(Wd) 


where, P s (x,y,zj) is the probability of survival through cell (x,y,z) in the yth direction to an 
adjacent cell. In the actual implementation, the values assigned to each cell are negative 
logarithms of the probabilities of survival. The problem is thereby transformed from one 



9 


of maximizing the product of survival probabilities to minimizing the sum of negative 
logarithmic probabilities. The probability is a function of terrain masking, fuel constraint, 
or time constraint. 

The investigators in artificial intelligence area [10] have suggested to use the heuristic 
search to find a near-optimal routes for autonomous helicopter. Two of the most 
commonly used heuristic search techniques for finding optimal path are the branch-and- 
bound and the A*-algorithm, discussed in Reference 10. The branch-and-bound is an 
exhaustive search method similar to both depth-first and breadth-first schemes. They 
search all possible paths until the goal is found. A*-algorithm is the branch-and-bound 
search in conjunction with the dynamic programming principle to reduce computations 
[28]. 

All these approaches employ the discretization of the terrain spatial coordinates before 
carrying out the search for the optimal trajectory. As a result, they assume that the route 
consists of straight line segments. On an uneven terrain, this implies that a large number of 
discretization intervals will be required to generate sufficiently smooth trajectories. 
Unfortunately, this increase in the number of discretization intervals is accompanied by an 
enormous increase in computational complexity. For example, in the case of discrete 
dynamic programming, this is of order {(n+l) 2 +l], where n is the number of 
discretization intervals in one spatial direction [12]. A solution advanced by some 
researchers for handling this "curse of dimensionality" is the use of parallel-computing 
architectures [30]. 

This report will propose alternative trajectory planning schemes based on the Euler- 
Lagrange equations [12]. These approaches require a one-dimensional search to determine 
optimal trajectories. Further details will be discussed in Chapter HI. 
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2.3 Previous Research on Helicopter Pursuit-Evasion 

Unlike the one-sided trajectory planning problem, reported research on two-sided 
trajectory planning has been very sparse. A previous research [31] employed a discrete 
matrix game approach for generating maneuvering decisions for low altitude flying 
helicopter during one-on-one air combat over a hilly terrain. Each player had seven 
maneuvering strategies, and thus the game matrix consisted of 49 payoff elements. Each 
element in this matrix represented the score evaluated using a scoring function. Under the 
perfect information assumption, the scoring function was composed of an orientation, a 
relative range, a velocity, and a terrain profile. The state variables required in evaluating 
the scoring function were obtained by numerically integrating the equations of motion for 
each of the seven strategies of the participants. After numerical integration, the saddle point 
was searched and optimal maneuvering strategies for each player were obtained. This 
procedure was repeated until terminal conditions are satisfied. 

According to Von Neumann and Mongenstem [32], every finite and discrete game can 
be cast in the matrix form. However, the dimensions of this matrix will be astronomical 
except for very simple problems. Additionally, the computational effort in conducting a 
search for the optimum can be prohibitive. R. Isaacs [17],provided the framework for 
obtaining solutions to continuous games with differential constraints. This will be further 
elaborated in Section 2.5. 

In this report, the helicopter pursuit-evasion problem will be studied as two one-sided 
optimal control problem using differential game theory. Two different formulations will be 
discussed. The first one requires two-dimensional search to determine optimal strategies. 
Another approach using nonlinear transformation techniques demands the specification of 
the terminal time. Details of these approaches will be given in Chapter IV. 



2.4 A Review of The Opt imal Control Theory 


In Older to motivate subsequent development, this section will present a review of the 
central results from optimal control theory [12]. 

Given the state equations: 

* = f(x,u,t) , x(to) = x 0 (2.5) 

where, x(t) : = state vector of dimension n, x e X 

u(t) : = control function of dimension m, u e U 

Initial constraints: 

R(x(t 0 ),t 0 ) = 0 

Terminal constraints: 

P (x(t,),t f ) = 0 and t f is free 


( 2 . 6 ) 


(2.7) 


Performance Index: 


J [u] = g(x(t f ),t f ) + J L(x,u,t) dt g) 

*0 

The optimal control problem is to pick u(t) to minimize J[u] while satisfying the state 
equations and the boundary constraints. 

The optimal control can be obtained using Dynamic Programming [27] or Pontryagin's 
Minimum Principle [33]. For most problems encountered in applications, these two 
approaches can be shown to be equivalent [12]. 


2.4.1 Continuous Dy namic Programming 


Define the continuous optimal return function as 

J°[x(tXt]= min J[x(tXu(tXt] < 2 - 9 ) 

ueU 

To simplify presentation, it is assumed here that the terminal cost is zero. Next, assuming 
the optimal return function to be continuous, one can write [34] 

t+e *f 

J°[x(t),t] = min { j L(x,u,x)dx + f L(x,u,x)dx } (2.10) 

u e U , J 

1 t+£ 

for sufficiently small £. If the vector functions u(x) and L are both continuous at t, there 
exist an e sufficiently small such that expression (2.10) can be approximated as: 

tf 

J°[x(tXt]= min {eL(x,u,x) + jL(x,u,x)dx) (2.11) 

ueU t+e 

From the definition of J° the optimal return function (2.9), this amounts to 

J°[x(t),t] = min {eL(x,u,t) + J°[x(t+E),t+e] } (2. 12) 

ue U 


The state evolution may next be approximated by 

x(t+e) = x(t) + ef(x,u,t) ( 2 . 13 ) 

Substituting equation (2.13) into (2.12), for sufficiently small positive e, and retaining only 
the first-order terms, one has 
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J°[x(t),t] = min {eL(x,u,t) + J°[x(t)+£f(x,u,t),t+e] } (2.14) 

U e u 

Taking a Taylor's series expansion of J [x(t)+ef(x,u,t)» t+E], and retaining only the first 
order terms, 

J°[x(t)+ef(x,u,t),t+e] = J°[x(t),t] + ef [x(t),t] f(x,u,t) + e J?[x(t),t] (2. 1 5) 

Next, substituting (2.15) into (2.14), and cancelling the J°[x(t),t] term and dividing by e, 
one has 

J°[x(t),t] = - min { L(x(t),u(t),t) + f[x(t),t] f(x,u,t)} (2.16) 

ue U 

Let u° (x,t) denote the optimal control given x and t. This control must yield a minimum 
for the right hand side of equation (2.16). Thus, 

T t [x,t] = - L(x,u°,t) - J° x [x,t] f(x,u°,t) (2.17) 

Equation (2.17) is known as the Hamilton- Jacobi-Bellman (HJB) equation. This is a first 
order nonlinear partial differential equation. 

This equation is difficult to solve if the functions L and f are highly nonlinear. As a 
result, this equation is often solved using the method of characteristics. The characteristics 
of the HJB equation are called the Euler-Lagrange (E-L) equations [12]. The E-L equations 
are first order nonlinear ordinary differential equations with prescribed boundary 
conditions. In the following we will indicate the derivation of E-L equations using the HJB 
equation. 
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Differentiating the expression (2.17) partially with respect to x, and noting that the 
control variables are independent of the state variables results in the expression 

£[x,t] = - L x (x,u,t) - J°„[x,t] f(x,u°,t) - f x (x,u ,t) r x [x,t] (2. 1 8) 


Now, the total derivative to J° x [x,t] is given by 


dJ*.[x,t] 

dt 


= J°, l [x ) t]+r iI [x,t]f(x,u,t) 


(2.19) 


Substituting (2.18) into (2.19), one has 
dJ° x [x,t] 


dt 


= - L x (x, u°, t) - f x (x, if, t) J° x [x,t] 


( 2 . 20 ) 


This set of first-order ordinary differential equations for J x °[x,t] can be solved if x and u° 
were known for all t and initial conditions for J x °[x,t] were given. J x °[x,t] are called the 
costates of the systems, often denoted by the variable X. 


2.4.2 Pontrvagin's Mi nimum Principle 

Introducing a new variable called the Hamiltonian [12], 

H (x,u,X,t) = L (x.u.t) + X T f (x,u,t) (2.21) 

The Euler-Lagrange equations (2.5), (2.20) can be written as 

x = f(x,u,t), ( 2.22a) 

X = -H x 


(2.22b) 


c 
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The fact that optimal control u has to minimize the quantities within braces in (2.16) leads 
to the so called optimality condition 

H u = 0 < 2 - 23 ) 

Equations (2.22) must satisfy the given boundary conditions. For additional details on this 
problem, see Reference 12. 

While the satisfaction of HJB is sufficient for optimality, additional conditions must be 
imposed while solving the optimal control problem using the E-L equations. In this case, 
the optimal controls emerging from (2.22), (2.23) should additionally satisfy the following 
condition [12]: 


(i) Legendre-Clebsch condition 

Huu> 0 

(ii) Weierstrass condition 
8H (x, X°, u°, u, t) > 0 


(2.24) 


(2.25) 


(iii) Jacobi condition 

Nonexistence of conjugate point in (tQ,t f ) 

Conditions (i) and (iii) are necessary for weak local minimum, while condition (ii) is 
necessary for strong local minimum. Strengthening conditions (i) and (ii) and closing the 
interval in condition (iii) constitute the sufficient condition [12]. In some situations, 
normality condition [12] has to be verified before testing the conditions (i) - (iii). 
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Additional conditions can be obtained using various combinations of these necessary 
conditions [35]. 


2.5 A Review of Differential Games 

If optimal control theory briefly reviewed in the previous section can be considered a 
theory for one-sided control problems, differential game theory may be identified as a 
theory for two-sided control problems. It has been shown [36] that the problem definitions 
and solution methods used in optimal control theory can be extended into the game theory. 

The theory of differential games is a subject concerned with the optimization of 
dynamic systems involving two or more players with conflicting interests. The study of 
differential games was initiated by Isaacs in 1954. In 1965, Isaacs published a book which 
details various aspects of differential games [17]. In 1957, Berkovitz and Fleming [37] 
solved a simple differential game using the calculus of variations. In a later research, 
Berkovitz treated a wider class of differential games using the calculus of variations [38]. 
Friedman's book [39] discusses the necessary conditions in differential games in terms of 
the more familiar optimal control theory notation. 

The aforementioned differential games mostly dealt with problems of the pursuit- 
evasion type having the zero-sum property, i.e., one player's loss being the other player's 
gain. Dropping the zero-sum hypothesis adds both conceptual and analytic complexity, but 
it may extend the utility of the theory of differential games to a much wider class of 
applications. Two typical non-zero-sum games are the Nash game and Stackelberg game, 
in which each participant has its own performance criterion. 
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Assuming that the players have perfect information of the current states and that their 
respective roles are determined before the game begins, a differential game may be stated as 
follows: 

Given differential constraints: 

x = f (x,4»,\|/,t), x (to) =x 0 

where, x(t) : = state vector of dimension n, x e X 

<]>(t) : = control of Player 1 of dimension i, <|> e ® 

\|r(t) : = control of Player 2 of dimension m, xjf € 'P 

Initial constraints: 

R(x(t o ),t o ) = 0 

and terminal constraints: 

P(x(t f ),t f ) = 0,t f isfree 

The terminal constraints define the stopping condition for the game. For example, if the 
participants have a "capture set", the game terminates at the instant the adversary enters the 

capture set. 

The performance index or payoff for the i th player may be defined as 

Ji[x,<|>,y,t] = gi(x(tf),t f )+ f Li(x,4>,Y,t) dt (2.29) 

Ju> 

Each player involved in the game attempts to optimize its own performance index with due 
attention to the other player's state-control variable evolution. Once the differential 
constraints, initial and terminal constraints and the performance index are defined, one may 
describe three different game scenarios. 


(2.26) 


(2.27) 


(2.28) 


2.5.1 Nash Non-Zero-Sum Game 


If an optimal solution exists, optimal closed loop control pair (<J>°, y°) must satisfy the 
Nash inequalities condition [40]: 




(2.30) 


J 2 (4»V) £J 2 (4>,V°). V <}»€<D 


(2.31) 


Inequalities (2.30) and (2.31) imply that optimal strategies for each player should yield the 
smallest cost for individual participants. Any deviation from the optimal strategy will yield 
a higher cost. 

2,5,2 Stacfaltorg Non-Zero-Sum gam? 


In this game, one assumes that the second player is the leader while the first player is 
the follower. As a result, the first player is operating a purely reactive fashion. If there 
exists a mapping M:*P — » <X>, and the following conditions are satisfied, then the pair (4>*, 
y*) e <J> x is called a Stakelberg strategy pair with Player 2 as a leader and Player 1 as 
follower [41]: 

Ji(My,y) < Ji(<t), v) ' (2.32) 

W,\i^)<J 2 (M W ) ' (2.33) 

<j>* = (2.34) 

In other words, the Stackelberg strategy is the optimal strategy for the leader when the 
follower reacts by playing optimally. An interesting property relating the Nash and 
Stakelberg strategies can be derived from equations (2.30) - (2.34) as follows [41]: 
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J 2 (<t>*,V*)^J2(<t> 0 .Y°) (2.35) 

which means that the leader in the Stackelberg solution achieves at least as good a cost 
function as the corresponding Nash solution. 

2.5.3 Zero-Sum Gam g 

Zero-sum games result when the two decision-makers are adversaries. One decision- 
maker's loss is the other decision-maker's gain. In this case, the equilibrium solution has 
the property that 

Ji =-J 2 = J (2.36) 

The above Nash inequalities criterion (2.30) and (2.31) may be reduced to 

J [x,«t>°,\|/,t] < J [x,<J>°,v°,t] < J [x,<|>y\t] (2.37) 

where, J [x,<t>°,v°,t] = V [x^y is called the value of the game. In such a differential game 
both players have the same performance index, with the first player minimizing it while the 
second player attempts to maximize. Equation (2.37) suggests that if the minimizing player 
deviates from his optimal strategy, the game will have a higher game value. Alternatively, 
if maximizing player deviates from the optimal strategy, the game will have a lower value 
than if the two were employing optimal strategies. This is the well-known saddle point 
condition in the game theory [32]. A similar saddle point condition may also be derived 
from the Stackelberg inequality conditions. 

Introducing the variational Hamiltonian 

H (x,<t>,y,X,t) = L + X 1 f 


(2.38) 
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and define the terminal conditions as: 

Q(x(t f ),t f ) = g(x(t f ),t f ) + v T P (xftj),^) 

Since * = f(x,u,t), one may write the performance index as 
J a = Q(x(t f ),t f ) + f (H - X T x) dt 

Ju> 


(2.39) 


(2.40) 


Since the control terms appear only in the function H, the saddle point condition (2.37) can 
be written as [42]. 


| H(x,<J>°, V ,X,t) dt <| H(x,4> c V,X,t) dt <| H(x,<|),v 0 ,X,t) dt 


(2.41) 


If for all t, functions L and f are continuous, the inequality (2.41) can be changed to 
the pointwise form [42] by using the principle of optimality which requires that at every 
instant controls should be chosen to make system optimal: 


H(x,4>°,v,X,t) < H(x,<J> 0 ,Y°,X,t) < H(x,<J>,\|/°,X,t) 


(2.42) 


Equation (2.42) is a sufficient condition for the inequality (2.41) and a differential game 
version of Weierstrass condition in the optimal control. 

Generally, by the order of action it is known [35] that 


max min H(x,<t>,y, X,t) = min max H(x,<l>,V,X,t) 
V + 4 ¥ 


(2.43) 
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Additional conditions sufficient for satisfying the above equation (2.43) as equality for a 
special case was given by Von Neumann and Morgenstem [32]. This is known in the 

literature as the Isaacs principle. 

If H is separable in <(> and y, i.e.. 


f(x, <J>, y, t) = fi(x, <(>, t) + f 2 (x, y, t) 
L(x,<j>,y,t) = Li(x> 0, t) + L 2 (x, y, t) 


(2.44) 

( 2 . 45 ) 


it may be shown [39] that 

max min H(x,0,y, X,t) = min max H(x,0,y ,X,t) = H(x,<{> ,y ,X., t) (2.46) 
V 4> <t> V 

Since the stationarity conditions (first-order necessary conditions) are the same for 
minimization or maximization, the differential game can be defined as a two-sided optimal 
control problem with coupling appearing via the transversality conditions. In this case, the 
necessary conditions for 4>(t) and y(t) to be optimal are: 


1) Euler-Lagrange equations 

X = f(x,<j>,y,t) (2.47) 

X = - H x (2.48) 

2) Transversality conditions at t f 

X T (t f )=Qx(x(t f ),t f ) (2-49) 

P (x(t f ),t f ) = 0 (2-50) 
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H(tf) = 0 , tf is free (2.51) 

3) Optimality conditions 

H 0 = 0 (2.52) 

H V = 0 (2.53) 

4) Isaacs principle 

max min H(x, X,t) = min max H(x, <t>,\|/, A,t) (2.54) 

VO 0 V 


It may be shown that this solution can be also obtained by solving a partial differential 
equation similar to the HJB equation [39]. In the differential game context, this PDE is 
called the Hamilton- Jacobi-Isaacs equation for the optimal J [39], viz, 

J°(x,t) = - min max H (x,<>,\j/; J°,t) (2.55) 

J°(x(t f ),t f ) = g(x(t f ),t f ) (2.56) 

The similarity of the differential game and optimal control is apparent from the 
foregoing. However, it is important to note that certain differences exist between optimal 
control problems and differential games [43]. First, although feedback control is not 
essential in the optimal control problems it becomes the central requirement in' differential 
games. Secondly, the solution is characterized by regions where solutions exist and where 
they do not. Moreover, verification of the second-order necessary conditions, while not 
routinely employed in optimal control, becomes mandatory in differential game solutions. 



2.6 Transformatio n of Nonlinear Systems 


Even though differential game theory has several potential applications in economics, 
military, and engineering, it has not been employed to the same degree as optimal control 
theory. However, linear-quadratic differential games have received considerable interest in 
the differential games literature [40, 41, 43, 44]. These games have been important in 
studying the local behavior of certain nonlinear differential games. Recently, it has been 
shown that a class of nonlinear differential games may be solved in closed form using 
transformations [18]. In that work it was shown that nonlinear transformation techniques 
are useful for implementing linear-quadratic differential game solutions in nonlinear 
differential games. 

2.6.1 Kronecker Indices and Brun ovskv Form 

Before tackling the nonlinear system control problem, a brief review of linear system 
theory in the state space will be given. For a given multivariable system 

x(t) = A x(t) + B u(t) , x e SR n , u e (2.57) 

the first step in designing control laws is to test for controllability. This involves the 
computation of a controllability matrix C[A, B} [45] as 

jbj, Ab v a\, .... A krl bj, b 2 , Ab 2 , ,A kafl b m J ' (2.58) 

For the system to be controllable, this matrix must have a rank n [45]. Next, this matrix 
may be normalized to determine a transformation matrix 


T = {en, ei 2 , .... eikp e 2 i, 


>®mli •• ®mk m ) 


(2.59) 
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Defining a new set of state variable x = Ty and substituting in (2.57), one has 


y = T- 1 ATy + 'T 1 Bu 


(2.60) 


or. 


y = Ac y + B c u 


(2.61) 


with 


A C = T 1 AT , B c =T 1 B 

Under this transformation the matrices A c and B c will have mostly zero or one entries 
together with a specific structure. For example, n=10, m=3, ki=3, k2=3, k3=4 




ri • •] 

1 000000000 


0 0 0 

01 00000000 


0 0 0 



0 1 • 

00 01 000000 

B„ = 

0 0 0 

00 001 00 00 0 

c 

0 0 0 



0 0 1 

00 00001000 


0 0 0 

00 00 0001 00 


0 0 0 

.00 00 0 00 01 0. 


.0 0 0. 


where symbol • means any numeric entry. In this controller canonical form, the canonical 
invariant parameters are 



(2.63) 


where, kj is called the controllability indices, or Kronecker indices and otj eigenvalues of 
system. 

Using control law 


u(t) = G v(t) - K y(t) 


(2.64) 


the system (2.61) becomes 
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y = (A e - B C K) y + B e G v 


(2.65) 


By suitable choice of input transformation G and state feedback gain K, the entries marked 
• in the {1st, (ki+l)st, (ki+k 2 +l)st,....} rows of Ac and B c can be zeroed out to get the 
special controller form as follows: 

y = A* y + B* v (2.66) 


For the aforementioned example (2.62), 


'0 

0 

0 

0 

0 

0 

0 

0 

0 

O' 


'I 

0 

0 ' 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

B* = 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 . 


0 

0 

0 . 


(2.67) 


As we can see in (2.67), all eigenvalues of system can be zeroed out. The only 
remaining invariants are the Kronecker indices. This special controller canonical form is 
called the Brunovsky canonical form [45] and exhibits a parallel array of m decoupled 
subsystems of dynamic order ki, k 2 , k 3 , k m [46]. 

2.6.2 Nonlinear Tra nsformation 

In the differential geometry, the Brunovsky canonical form can be viewed as a group 
acting on the space of linear systems. Each subsystem forms a single orbit under the group 
action and new system pair {A*, B*} are cyclic (single-companion-matrix). Such 
differential geometric concepts in linear control theory have been extended into nonlinear 
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systems with control variables appearing linearly in the dynamics by Krener [47] and 
Brockett [48], An infinitesimal group theory and tangential transformations for nonlinear 
differential equations studied by Marius Sophus Lie [49] forms the basis for transforming a 
nonlinear system to linear system. Compared with the usual linearization using a Taylor 
series expansion about a fixed equilibrium point, local tangential transformations expand 
Taylor series continuously along the trajectory. Thus, the approximation, like a turtle, 
carries its convergence house with it [50]. 

Consider the nonlinear system 

m 

^ = fOO ^ gi(x)Ui(t) ^2 ggj 

where x, f(x), g,(x) are n vectors with the hypothesis that f(0) = 0, causality condition, and 
u are the m control variables. This system may be transformed to Brunovsky's canonical 
form using two distinct approaches. Each of these are discussed below. 

2.6.2. 1 Hunt & Su's Approach 

References 51 and 52 showed that the transformation T = (Tj, T 2 ,...., T n , T n+1 ,...., 
Tn+m) is required to have the following properties 

(i) T (0) = 0. 

(ii) Tj, T 2 ,...., T n are functions of x x , x 2 ,...., x n . 

(iii) T n+1 , T n+2 T n+m , are functions of x lf x 2 ,...., x^Uj, u 2 , u m . 

(iv) T maps the open set U of 9t n into 9t n with a nonsingular Jacobian matrix. 

(v) Tj, T 2 ,..„, T n are the state variables and T n+1 , T n+2 , ...., T n+m are the controls 
for a linear time-invariant system in the Brunovsky form. 

Reference 53 gives necessary and sufficient conditions to map to the Brunovsky form 
with m Kronecker indices Kj, k 2 , k^. Because the Lie bracket operation on pairs of 
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vectors keeps invariant characteristics independent of the choice of coordinate systems 
used, introduce the Lie brackets 

< 2 - 69 > 


where 3g/9x and df/dx are n x n Jacobian matrix. One may define an alternative notation to 
simplify the analysis. Thus, let 


[f,g] = (ad'f.g) 


(2.70) 


Successive Lie Brackets can then be expressed as 
g = (ad°f,g) 

[f,g] = (ad^.g) 

[f,[f,g]]= (ad 2 f,g) 

[f,(ad 2 f,g)] = (ad 3 f,g) (2.71) 

[f,(ad nl f,g)] = (ad^.g) 

The transformation T = (Tj, T 2 ,...., T n , T n+1 ,...., T n+m ) exists if and only if 

(i) the set C = {(ad°f,gi), (ad^.gi) (ad Kl ‘ 1 f,gi), (ad°f,g 2 ), (ad 1 f,g 2 ) 

(ad K2 l f,g 2 ), ..... (ad°f,g m ), (ad^.gm), (ad^f.gn,)} spans R n about the origin. 

(ii) the sets Cj = {(ad°f,gi), (ad'f.gi) (ad K i' 2 f,gi), (ad°f,g 2 ), (ad 1 f,g 2 ),...., 

(ad K J' 2 f,g 2 ), (ad°f,g m ), ..., (ad K f 2 f,g m )) are involute for j = 1, 2 m. 


(iii) the span of each Cj is equal to the span of Cj o C. 
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By definition, a linearly independent set of vector fields {x lt xj is involute if and 
only if there are scalar functions such that 



for all i,j,k 


(2.72) 


2.6. 2.2 Mever's Approach 

This approach is more intuitive. In this technique, the states to be controlled are 
successively differentiated until control terms appear in the equations. Various steps during 
these differentiations form the mapping. 

2.6.2.3 Comparison of Two Approaches 

Although Hunt and Su's approach is more systematic than Meyer's, it is nearly 
impossible to solve partial differential equations for T-transformations except in very 
simple problems. Meyer's approach is a special case of Hunt and Su's approach. 

If system equations are derived by the classical dynamics with forces and moments as 
control variables, then control terms will appear in the second-order kinematic equations. 
In this case, Meyer's approach can be easily implemented and the transformed system is 
always a double integrator system. This technique has been used in Robotics for several 
years and called the Computed Torque Method. 

For comparisons, a simple example is provided by the problem of a vertically moving 
vehicle of unit mass which has thrust u and drag being proportional to the square of the 
speed. The equation of motion can be written 

x = -G-Kx 2 +u (2.73) 


with G denoting gravity. In state variable form, the system can be expressed as follows: 
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*i = x 2 

X 2 = - G - Kx 2 + u 


In the standard notation (2.68), it can be expressed as follows: 


f = 


X2 

-G-K* 



(2.74) 

(2.75) 


(2.76) 


First, consider Hunt and Su's Method. The necessary and sufficient conditions are 
that rank {g, [f,g]} =2 and that {g, [f,g] } be involute. It is easy to check that 


[f,g] = 0 - 


0 

0 



- 2Kx 2 


(2.77) 


g. [f.g]| = 


0 

1 


- 1 

2Kx 2 


1*0 


(2.78) 


which has rank 2 and the vector field {g, adf(g)} is involute. Then, the method begins 
with writing a Brunovsky form 



(2.79) 


where 


yi = Ti(xi,x 2 ) 

(2.80) 

yi = T 2 (xi,x 2 ) 

(2.81) 

v = T3 (xi,x 2 ,u) 

(2.82) 


From Brunovsky form. 
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?i = ti(xi,x 2 ) = *1 + * 2 = <dTi,f> + <dTi,g>u = T 2 (xi,x 2 ) (2.83) 

y 2 = t 2 (xi,x 2 ) = |^-*i + |^-* 2 = <dT 2 ,f> + <dT 2 ,g>u = T 3 (x lt x 2 ,u) (2.84) 


with 


<dT, f> = fi + f 2 
dxj 0 X 2 


(2.85) 


From (2.83) and (2.84), one has 

<dTi, f> = T 2 , <dTi, g> = 0 (2.86) 

<dT 2 ,g> * 0 (2.87) 


Using Frobenius Theorem [54], i.e., 

<dh,[f,g]> = <d<dh,g>,f> - <d<dh,f>,g> (2.88) 


one can get a following relation 

<dT lt [f,g]> = <d<dT 1)g >,f> - <d<dT!,f>,g> = 0 - <dT 2 ,g> * 0 (2.89) 


Therefore, from (2.86) and (2.89) 

0X2 


<dT lt [f,g]>= .*±*0 


The simplest solution is 


(2.90) 

(2.91) 


yj =Ti =Xi 


(2.92) 



and from (2.86) 


3 1 


y 2 = T 2 = <dTi,f> = x 2 


(2.93) 


The feedback linearization control u can be obtained from (2.84) as follows: 

v - <dT 2 ,f> s v + GjJCx^ = y + G + 

<dT 2 ,g> 1 


(2.94) 


Second, consider Meyer's approach. Since the motion is single degree of motion, 
change coordinate by setting 


Y = x 


(2.95) 


After differentiating until control term u appears, the system can be expressed as follows: 


Y= x 


Y = X 2 = -G-Kx^ + u 


(2.96) 


(2.97) 


Define 


Y = v 


(2.98) 


The Brunovsky controller form is 

*ii =r° 1 

yj Lo o 


>i 

y % j 


(2.99) 


with the feedback linearization control 


u = v + G + Kx^ 


( 2 . 100 ) 
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2. 6. 2.4 Applications 

References 55, 56, and 57 designed an automatic flight controller for UH-1H 
helicopter applying nonlinear transformation techniques and showed that the controller 
exhibits good performance in all flight modes. As discussed in previous section, nonlinear 
transformation technique based on Meyer approach needs to successively differentiate the 
controlled states to get transformation map from nonlinear system to linear system. Meyer 
presented a numerical approach for calculating these derivatives. A successive numerical 
differentiation, however, requires formidable calculation, with attending numerical 
difficulties. 

To reduce the amount of computations, Menon [58] introduced singular perturbation 
technique to simplify the nonlinear mapping for a flight test trajectory controller of high 
performance airplane. The slow time scale controller follows path command and generates 
steady state values for the body attitude. The fast time scale controller is designed to track 
the commanded values for the body attitude. In Reference 18, the nonlinear 
transformations have been used to derive pursuit-evasion guidance laws for high 
performance aircraft 

Heiges [59] applied the above mentioned techniques to design a helicopter trajectory 

* 

controller and implemented it on the TMAN program. This controller has demonstrated 
excellent performance in various tactical flight modes. 



CHAPTER m 


TRAJECTORY PLANNING AS A ONE-SIDED 
OPTIMAL CONTROL PROBLEM 


3.1 Introduction 

This chapter discusses the problem formulation and optimal trajectory synthesis using 
two diff erent performance indices. Candidate trajectories are generated and their optimality 
is tested using second-variation analysis. Numerical effort involved in the computations 
are analyzed. 

For the first Optimal Route Planning problem, the terrain constraint is embedded into 
state equations via a coordinate transformation. The performance index is a linear 
combination of flight time and terrain altitude. In this problem, the computational algorithm 
requires the first and second partial derivatives of the terrain profile. 

The second Optimal Route Planning problem uses a performance index consisting of a 
quadratic form in the terrain altitude, lateral deviation from the nominal trajectory, and 
heading angle. By changing the independent variable from time to down-range, the order 
of the problem is reduced. Each of these trajectory planning schemes will be discussed in 
greater detail in the ensuing. 


A kinematic model of the helicopter will be employed in the ensuing analysis. Let the 


terrain profile be specified by a function 

h t = f(x,y) (3.1) 

where h t is the altitude above a preselected datum at any specified position (x,y), x and y 
being the down-range and cross-range measured in a chosen inertial frame. It is assumed 
here that h t > 0 and that the terrain f(x,y) has continuous first and second partial 
derivatives. This fact is important to ensure that the trajectories emerging from this optimal 
trajectory planning problem are implementable. Additionally, this is consistent with the 
proposed cubic spline parameterization of the digital teiTain data. While executing nap-of- 
the-earth flight, the helicopter altitude motion will follow the terrain profile (3.1) with a 
specified altitude clearance. As a result, the helicopter altitude at any location (x,y) is given 
by the equation 

h = f (x,y) + h c (3.2) 

In (3.2), h is the helicopter altitude and he is the specified terrain clearance. For NOE 
flight, the clearance is between 5 and 120 feet [1]. 

A sample terrain profile with the x, y, h coordinate system is shown in Figure 3.1. 
The local coordinate system x € , y e , z t used in subsequent analysis is also defined in this 

figure. This moving coordinate system has its origin on the terrain surface at the current x, 
y position with x^-y^ plane being the tangent plane. The principal direction of this system 
is along the intersection of the x^-y^ plane with the x-h plane. Accordingly, z ( points in the 
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direction of the normal vector to the local tangent plane. The transformation of vectorial 
quantities from one system to the other can be accomplished using the terrain gradients: see 
Appendix A for details. Since the helicopter is constrained to move on the surface defined 
by equation (3.2), the velocity vector lies in the instantaneous x € -y^ plane, with the angle % 

defining its orientation on this plane. The helicopter velocity components in the local frame 
can be resolved as: 


a* = v cos x 

(3.3) 

Si = V sin X 

(3.4) 


The local heading angle x and the airspeed V are assumed to be the control variables in the 
present trajectory planning problem. 

Note that it is important to include velocity as a control variable in the present problem 
to ensure hodograph convexity required for the existence of optimal controls [60]. In order 
to ensure that the control emerging from this formulation are implementable, the helicopter 
speed is next bounded as 

0<V mm < V<V max (3.5) 

Because a simple kinematic system is considered here, V mW i corresponds to the speed at 
which sufficient excess power is available for avoiding unknown obstacles. 

The velocity components (3.3) and (3.4) may next be transformed to the down-range, 
cross-range, altitude frame using the relations 
_ V COS X | V f*fy sin X 

Vl+fx V l+f*V 1+ft+fy 


(3.6) 
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- V V 1+f \ sin x 
V 1+fi+fy 


(3.7) 


The altitude rate on the terrain is given by: 

h = f x x + f y ? (3.8) 

These relations are derived in Appendix A. In the expressions (3.6)-(3.8), f x and f y are 
local terrain gradients, assumed to be calculated from the given terrain profile (3.1). 

In equations (3.6) and (3.7), the quantity 1/V 1+fi+fy denotes the cosine of the angle 
between the vertical z-axis in the inertial coordinate system and the direction of the normal 
in the local tangent plane. And the quantity 1/Vl+ff denotes the cosine of angle between 
the down-range x-axis in the inertial coordinate system and the x t axis in the local tangent 

plane. 

It is sometimes desirable to include ambient winds in the trajectory planning problem. 
If the winds aloft along down-range and cross-range directions are given by 

u = Q(x,y), v = R(x,y) (3.9) 


these may be added to the right-hand side of (3.6) and (3.7) to obtain the equations of 
motion as 


V cosx 

VI*8 


+ 


V f x f y sin % 
Vi+f£V l+fi+fy 


(3.10) 


- V Vl+ff sin X | v 
V 1+fi+fy 


(3.11) 



e 
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The equations of motion (3.8), (3.10), (3.11) may be used whenever the helicopter is 
flying in a terrain-following/terrain-avoidance mode. 

Known threats and obstacles may be incorporated in the trajectory planning problem 
by defining threat overlays of the form 

Ah = P(x,y) (3-12) 


and adding them to the basic terrain profile given by equation (3.1). The composite profile 
may then be used to define the equations of motion (3.10) and (3.11). In that case, the 
resulting trajectories will exhibit automatic threat and obstacle avoidance characteristics. 

Additionally, it is possible to consider a formulation in which the specific energy of the 
helicopter is maintained constant. This will occur whenever the throttle is set to maintain 
thrust equal to drag while executing the nap-of-the-earth flight. In this case, the airspeed 
will depend on the terrain profile as: 


V = Y2g [E - f(x,y) - he] 


(3.13) 


In (3.13), g is the acceleration due to gravity and E = h + V 2 /2g is the specific energy. 



The performance index considered in this problem is a linear combination of flight time 
and a terrain masking function. Following the existing literature [7], trajectory masking 
will be assumed to be accomplished if an integral proportional to the helicopter altitude is 
minimized. Admittedly, this masking function is crude since it is based on the contention 
that depressed terrain tends to provide a better masking. If improved terrain masking 
functions given as a function of down-range x and cross-range y were available, they can 
be included in the following analysis without difficulty. For simplicity, in all that follows, 
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the teirain masking will be assumed to be accomplished if the integral of helicopter altitude 
is minimized. A relative weighting factor is next introduced between the flight time and the 
terrain masking function to control the trade-off between these two, often conflicting 
requirements. Thus, a composite performance index of the form 

J = J [ (1-K) + K f (x,y) ] dt (3.14) 

o 


with 


0 £ K < 1 


(3.15) 


will be used in the following. 
The initial conditions 


x(to) = x o- y( t o) = y<) : p ycn 


(3.16) 


and the terminal conditions 


x(tf) = x f , y(t f ) = y f : given 


(3.17) 


The variational Hamiltonian [12] may next be formed -by adjoining the differential 
constraints (3.6), (3.7) to the performance index (3.14) to yield: 

H=l-K + Kf(x,y) + X,M+ , Yl -^ ) ' 

Vi+f! Vl+flVi+fi+f? 


, ,-VVl+fl sin X s 

M Vi+fS+f? } 


(3.18) 


Note that the wind components have been dropped from (3.18). 
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The Euler-Lagrange equations for this optimal control problem are: 

B! sin % ky + (B 2 sin % + B 3 cos %) 


K = -Kf x - 


A\Al 


Xy = -Kfy 


B 4 sin x + (B 5 sin % + B 6 cos %) y 


a\a\ 


where 

Ai=Vl+f| 

A2 = Vl+ f| + fy 

Bl = {- A| fxfxx + A I (fyfxy + fx^xx)} Aj 

B2 = - [ fxfxx A2 + Aj(fyf X y + fxfxx) } fxfy + Aj A2 (f x fxy + fy^xx) 

B3 = - Aj fxfxx 

B4 = {■ A2 f x fxy + Aj (f x fxy + fyfyy) ) Aj 

B5 = - { f X f X yA^ + Aj(f x f X y + fyfyy) )fxfy + Aj A| (fyfxy + fx^yy) 
B 6 = - A\ fxfxy 


with the optimality condition 


tan x = 


Xxfxfy - >. y (1+g) 

X x V l+ft+f2 


(3.19) 

(3.20) 


(3.21) 


The equations (3.6), (3.7), (3.19), and (3.20) together with (3.21) constitute a nonlinear 
two point boundary value problem, which can be solved if the initial conditions on the two 
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costates X x and X y were known. However, since the variational Hamiltonian is not 
explicitly dependent on time and the final time is free, this optimal control problem has a 
constant of motion, viz., 

H(t) = 0, 0<t<t f (3.22) 


i.e.. 


0.1-K + Kf(x.y) + X,(^, Vf f in * ) 

Vl+f \ Vl+f?Vl+ft+fy 


. , /- V VI+f£ sin x. 


(3.23) 


This constant of motion may be employed to eliminate one of the costates in the 
problem. Using (3.21) and (3.23), one may solve for the costates Xx and Xy as 

K -Il-K + Kf(x,y))/T7f3 _ c _ (3.24) 

, { 1 - K + K f(x,y))(V 1 + f? + sin % - fxfy cos X) (3.25) 

A _ 

vVi + ft 


Additionally, the costates can be completely eliminated from this problem by employing an 
adjoint-control transformation as illustrated in the following. 

The expressions (3.24) and (3.25) are next differentiated with respect to time and 
equated to the right hand sides of the equations (3.19) or (3.20). This process yields a 
differential equation for x as: 

• {(A 1 K+A 2 )cosx+A 3 (A 4 K+A 5 )sinx} V 
X = 


A 6 (A 7 K + 1) 


(3.26) 



where. 


B x = (f-l)fxx - d+fx 2 ), 
B 2 = (l-f)(l+fx 2 ), 


B 3 = fx(l+fx 2 ) 2 , 


Aj — BjA3fy, 
A 2 = fxxfyA>3* 


a 3 =Vi + f£ + f| , 

A4 = Bjf x fy^ + B 2 f X yfy " B 3 , 
A5 = fx^xx^y 2 * (l+fx 2 )fxyfy, 

A 6 = {(l+f x 2 )(l+fx 2 +fy 2 )} 3/2 . 


A 7 = (f-1). 

An implicit assumption made in deriving (3.26) is that X exists everywhere on the 
terrain. This aspect will be verified while discussing the second variation analysis for this 
problem. The expression (3.26) was obtained using the MACSYMA program [61]. 

With the foregoing analysis, the optimal route planning problem has been reduced to 
that of solving a set of three nonlinear differential equations (3.6), (3.7), and (3.26) with 
one unknown boundary condition x(0)- The solution of this problem requires the 
determination of the initial value of heading angle X- Since x(0) and y(0) are known, X(^) 
must be selected such that the final conditions on x and y are the desired values x(tf) and 
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y(t f ). A simple iterative technique such as the method of bisections [62] may be set up to 
solve this problem. Flow chart of such an iterative scheme is illustrated in Figure 3.2. If a 
solution for the system (3.6), (3.7) and (3.26) satisfying the given boundary conditions 
exists within the given x(0) range, then it can be shown that the scheme given in Figure 2 
will find it in finite number of iterations [62]. Moreover, enforcing the conditions for the 
existence of optimal controls can yield further guarantees on the convergence of this 
numerical algorithm. 

Consider next, the second control variable in this problem, viz., the helicopter speed. 
Since the second control variable V appears linearly in the variational Hamiltonian and is 
bounded, the optimal control is given by 


V = V ma *, if S < 0 


V = V, 


nun* 


if S >0 


V = Singular, S ■ 0 


(3.27) 


where S is the switching function obtained from 


S = 


9H 

av 


namely, 

s _ K ( Vl + f£ + fj? cos % + fxfy sin x) - Xy(l + f£) sin % 

Vl +f? Vl + f \ + f£ 


(3.28) 


(3.29) 


Substituting X x andXy from (3.23) and (3.24) into (3.29) 
s _ [(1 - K) + K f} 


V 


(3.30) 


Since V is always positive, the sign of the switch function is determined by the term within 
the braces. This term is always less than zero by definition (0 < K < 1, 0 < f). This 
expression suggests that the maximum speed setting is optimal throughout the trajectory. 

Euler solutions for the optimal trajectory planning problem may be generated by 
numerically integrating the three first order nonlinear differential equations (3.6), (3.7), and 
(3.26). Starting from arbitrary initial conditions x(0) and y(0), Euler solutions to various 
end conditions can be generated by changing the initial value of the heading angle. In the 
present work, a sample terrain data from the U.S. Geological Survey [63] was used. The 
terrain approximates a part of the Nassau Valley area in California shown in Figure 3.3. 
The terrain data is stored at 1000 feet intervals and interpolated using Cubic Spline Lattices 
[64]. This terrain data is given in Table 3.1. First and second gradients of the terrain 
profile required in subsequent calculations are generated by differentiating the spline 
polynominals analytically and substituting for down-range and cross-range values. The 
nonlinear differential equations are integrated using a fixed-step fifth-order Runge-Kutta- 
Merson technique and the method of bisections is used to carry out the one-dimensional 
search. All computations were carried out on a VAX 11/750 computer system with double 
precision arithmetic. 

Figure 3.4 illustrates time-optimal trajectories starting at the point O and terminating at 
several end points. These trajectories were obtained by setting K= 0 and varying the initial 
value of the heading angle. This value of K corresponds to the case of time-optimal 
control. The trajectories appear to be close to straight lines except in regions of large terrain 
curvature. A family of Euler solutions starting at the point O with a large weight on the 
terrain masking (K = 0.99) is given in Figure 3.5. These trajectories exhibit a more 
significant curvature. An interesting feature of this solution family is that some of the 
trajectories appear to intersect in certain regions of the given terrain. This implies the 
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existence of more than one trajectory satisfying the stationary and boundary conditions. In 
this situation, the selection of a particular path has to be based on second order necessary 
conditions. Such an analysis will be presented in the Section 3.4. For a typical set of 
boundary conditions. Figure 3.6 illustrates the difference between time optimal and 
maximum terrain masking trajectories with initial heading angles 50 and 68 degrees, 
respectively. It may be observed from this figure that the temtin masking trajectory tends 
to seek out lower elevations while time optimal trajectories appear to minimize the arc 
length. 


3.3 Optimal Route Planning Problem No.2 

In the trajectory planning scheme discussed in the foregoing, the computation of 
extremals required first and second partial derivatives of the terrain profile. In some 
situations, it may not be desirable to compute these derivatives due to the nature of terrain 
data. A formulation of the trajectory problem that does not require these partial derivatives 
will be discussed next. 

3.3.1 Problem Formulation 

Assuming that the helicopter has a speed of V, the flight path angle y and the heading 
angle y, the velocity components in the defined inertial reference frame is given, by 

* = V cosy cosy (3.31) 

y = V cosy siny ( 3 . 32 ) 

The heading angle is the control variable in the this problem, while the flight path angle is 
defined by the terrain profile. This is due to the fact that the helicopter is executing terrain 
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following flight. Thus, the kinematic model of the helicopter flight is given by the 
differential equations (3.31), (3.32) and a nonlinear algebraic equation (3.2). In addition to 
this, one can define three differential equations describing the point mass dynamics of the 
helicopter. While it is desirable to include this in the formulation, the resulting trajectory 
optimization problem becomes intractable. Note that it is possible to correct the present 
results for neglected dynamics using singular perturbation theory [65 - 67]. The present 
research will not address this aspect of the optimal trajectory synthesis problem. 

In the ensuing formulation, time is not included in the performance index. Moreover, 
since time does not appear explicitly on the right-hand side of equations (3.31) and (3.32), 
the independent variable is next changed from time to down-range. This yields a dynamic 
equation of the form 

y' = tan \j/ (3.33) 

Here, a prime over the variables represents differentiation with respect to down-range, the 
independent variable in this problem. Note that this formulation is independent of the 
vehicle velocity V. As a consequence, it is possible to impose an additional acceleration 
constraint on the problem. It needs to be underscored that the vehicle velocity cannot be 
permitted to be zero along the trajectory. Otherwise, the present modelling will lose its 
validity. 

3.3.2 Optim al Trajectory Planning 

Assuming that the nominal trajectory to be flown by the helicopter is given by the 
function y c (x), the objective of the second trajectory planning scheme is to maximize terrain 
masking while minimizing deviations from the nominal trajectory. With this point of view, 
the equation (3.33) may be modified as: 
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6y' s tan y - y c '(x) (3.34) 

Here, y c '(x) is the derivative of the command y c (x) with respect to down-range x. The 
heading angle y is the control variable in this problem. The second optimal control 

problem is then defined as: 

x f 

min 4[ (f 2 + e5y 2 + ay 2 )dx (3.35) 

v(t) 

subject to the differential constraint (3.34). The quantities e and a are factors that control 
the relative weight between deviations from the specified path and lateral acceleration. For 
mathematical convenience, we next replace the term corresponding to y 2 with tan 2 y. 
Moreover, the nominal trajectory is often specified by straight line segments. In this case, 
one can redefine the origin of the coordinate system at the initial point with the abscissa 
pointing in the direction of the down-range direction without any loss of generality. In this 
case, the quantity 5y can be set equal to y. 

With these modifications, the optimal control problem is redefined as 

min 
y’ 

The Euler's necessary condition for this problem can be obtained [12]: 

_ (ey 1 ffy) 

y a (3.37) 

Equation (3.37) is a nonlinear second order differential equation with a varying coefficient. 
The initial condition y(x 0 ) and the final condition y(x f ) are specified. As in the previous 

section, the quantity f y is the gradient of the terrain profile in the cross-range direction. The 
Euler's necessary conditions can be obtained via two distinct, although equivalent 


-±J '(f 2 + ey 2 + ay' 2 )dx 

^ r\ 


(3.36) 
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approaches. First, one may use the classical calculus of variations to derive (3.37). 
Alternative approach is to let y" = U and proceed via modem optimal control theory [12]. 
Here, U is the control variable. 

Numerical solutions to the differential equation (3.37) yield the extremals. To 
construct an extremal joining a pair of x,y boundary conditions, the unknown initial 
condition y'(0) needs to be determined. 

Since just one unknown parameter is involved, it is possible to employ the method of 
bisections to find the solution. Moreover, linear interpolation may be employed for terrain 
interpolation since the method outlined here requires just the first gradient of the terrain 
profile. 

Figure 3.14 shows a typical set of Euler solutions starting at the point O for a = 10 5 , e 
= 0.001. These trajectories were generated by changing the initial value of y' and 
integrating the Euler's necessary condition forward. Note that the effect of increasing a is 
to produce trajectories that are closer to straight lines, while the influence of increasing £ is 

to introduce more features into the solution. 

Figure 3.15 illustrates two trajectories with same boundary conditions, one being a 

straight line joining the initial and final points and another being the extremal generated by 

integrating equation (3.37). This trajectory corresponds to a = 10 6 , £ = 0.01, and initial 

* 

heading angle = -45 deg. Nature of the extremal agrees with general intuition that 
helicopter should fly at lower altitudes and its path may zigzag around hills. Altitude 
profiles along two trajectories in Figure 3.16 further illustrate this fact 
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2A Second-Variation Analysis 

It can be shown that a sufficient condition for the extremals obtained in the foregoing 
to provide a weak local minimum for an optimal control problem is that the second variation 
be strongly positive [68]. The second variation will be strongly positive if the Legendre- 
Clebsch necessary condition is met in the strengthened form and the no-conjugate point test 
is satisfied. In optimal control problems with unbounded controls, it can be shown that the 
variational Hamiltonian is equivalent to the Weierstrass excess function [68]. As a result if 
the second variation tests are satisfied, then the extremals provide a strong local minimum. 
In addition to this, it is known [69] that if the integrand of the performance index satisfies a 
convexity condition, the extremals affording strong local minimum also provide a global 
minimum. In the following, each of these conditions will be examined for the two 
previous problems. 

3.4.1 Second-Variation Analysis for O.R.P. #1 


3.4. 1.1 Legendre-Clebsch Necessary Condition 

In an earlier section, it was shown that the optimal value of airspeed V = V miK . Since 
this control variable appears linearly in the Hamiltonian, theXegendre-Clebsch condition 
reduces to the scalar form 


h»20 


(3.38) 


Taking second partial derivative of the Hamiltonian (3.18) with respect to control %, 
-MVi + fj + f? cosx + f x f y sin%) + Xy(l + f%) sinx 


Hxx = ' 


Vl +fx Vl +f| + fy 


(3.39) 


i 
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Next, substituting X* and Xy from (3.24) and (3.25) into (3.39), one obtains 

H xx = (l-K) + Kf(x,y) (3.40) 

The right-hand side of (3.40) is strictly positive by definition. Thus, the Legendre- 
Clebsch condition is satisfied in the strengthened form at every point on the terrain. The 
implication of this is that the problem is regular. As a result, the extremals will be smooth 
and provide a weak local minimum for sufficiently short intervals. 


3.4. 1.2 Weierstrass Test 


In the Calculus of Variations, strong variations only bound the magnitude of 5x, while 
weak variations bound both the magnitude of 8x and the magnitude of the time derivative of 
5x [35]. Weierstrass gave a test for verifying strong local minimum characteristics of 
extremals. It can be shown that the Weierstrass's excess function [68] is equivalent to the 
Variational Hamiltonian in optimal control theory [12]. 


H(x,X,*,t) = l-K+Kf+X x V{ 


cosx , f x f y sin X ] 

Vl+ff Vl+ffV l+f x +fy 


X , V l+fjsinx 

Y VwN? 


Substituting X x and Xy from (3.24) and (3.25) into (3.41), one has 

H(x,X,x,0 = (1 - K + K f) { 1- cos (x° - X)) 

Here, x° is the value of x satisfying the optimal condition H^=0. 
and 


(3.41) 


(3.42) 


H(x,X,x,t) > H(x,X,x°,t) = 0 


(3.43) 
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In the foregoing it has been shown that the Weierstrass necessary condition is satisfied 
along the extremal. Thus, the performance along the extremal is lower than along any other 
trajectory. 

3.4.L3 Jacobi Test 

For extremals of finite length, the task of ensuring that the second variation is 
nonnegative for admissible variations leads to the accessory-minimum problem in the 
calculus of variations. This problem attempts to produce a system of admissible variations, 
not identically zero, which offer the most severe competition in the sense of minimizing the 
second variation. If a system of nonzero variations making the second variation zero can 
be found, then a neighboring trajectory is competitive. In this case, the test extremal 
furnishes at best an improper minimum and at worst a merely stationary value [70]. First 
value of the independent variable for which such a nontrivial system of variations can be 
found defines a conjugate-point 

It has been shown in die references [12] [71] [72] that the accessory minimum problem 
leads to the analysis of the nature of solutions to linearized Euler-Lagrange equations. We 
note here that this problem may be cast in the standard linear-quadratic format using the 
backward sweep method [12]. However, algebraic linearization of the Euler-Lagrange 
equations (3.6), (3.7), (3.19), (3.20) can become excessively involved. As an alternative, 
a numerical conjugate point test will be employed in the present research. Reference 73 
discusses several numerical methods available for conjugate-point testing. A direct 
approach for conjugate point testing will be pursued in the following. This approach is 
based on Theorem 12.1 -12.3 in Reference 70. This theorem is discussed in detail in 
Appendix B. 



From the existance theory of linear ordinary differential equations [74], it is known 
that the solutions to the linearized Euler-Lagrange equations should be linearly independent 
More concisely, the Wronskian determinant should not be zero on an open interval. For 
these fundamental solutions to be unique, they should be nontrivial and satisfy initial 
conditions. 

In the present problem, the fundamental solutions are obtained by using the finite- 
difference approximation 


8xi (t) = x*(t) - x^t), 5yi (t) = y*(t) - y*(t) 

dx 2 (t) = x*(t) - xi(t), dy 2 (t) = y*(t) - yi(t) 


(3.44) 


where (x^O.y'ft)) and (xi(t).yi(t)) are neighboring extremals generated by perturbing the 
initial value of heading angle. Note that this is equivalent to perturbing the initial value of 
costates, while maintaining the appropriate state initial conditions. 

The characteristic determinant A(t) is then formed as: 


8xi(t) 8x 2 (t) 
8yi(t) Sy 2 (t) 


(3.45) 


From the theory of differential equations [74] it is known that this determinant cannot be 
identically zero. If this characteristic determinant (3.45) after being zero at initial t = 0, 
subsequently becomes zero at t = t*, with t* < tj, then the point t* is said to be conjugate to 

the initial point t = 0. 

This numerical conjugate point test is applied to the maximum masking extremals 
given in Figure 3.5 with K=0.99. In this figure, the extremals A and B are of particular 
interest since these are competing extremals satisfying the boundary condition pair (O, F). 
The characteristic determinant (3.45) evaluated along these extremals is given in Figures 
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3.8 and 3.9. From these figures it is clear that a conjugate point to point 0 occurs along the 
extremal B, while none is encountered along extremal A. Thus, the extremal A affords a 
strong local minimum if the desired end conditions were points O and F. On the other 
hand, the extremal B provides merely a stationary value. This fact is confirmed by 
computing the performance index along these trajectories and given in Figure 3.10. This 
figure indirectly proves the princple of optimality by showing that the performance index 
along optimal trajectory A between two points O and F is always lower than neithboring 
extremal B. The point conjugate to the point O along extremal B is denoted by the point C 
in Figure 3.5. A rule of thumb given in Reference 75 is that a conjugate point can be 
expected to exist at point where the tangent to the extremal experiences a sudden change. 
In the present case, it may be observed that the extremal B experiences a sharp turn at the 
point C. 

In the minimum flight time problem, K = 0, the integrand of the performance index is 
a constant, and it does satisfy the weak convex condition given in Reference 69. Thus the 
extremals satisfying three second variation tests also provide a global optimum [69]. 

In Figure 3.11, two extremals satisfying the boundary condition pair (O, F) are given 
for the min-time criterion. Following the same procedure as for the maximum masking 
problem, the characteristic determinant evaluated along these extremals is given in Figures 
3.12 and 3.13. From these figures it is clear that a conjugate point to point O occurs along 
the extremal A, while none is encountered along extremal B. Thus, the extremal B affords 
a global minimum if the desired end conditions were points O and F. On the other hand, 
the extremal A provides merely a stationary value. 
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3.4.2 Second-Variation An alysis for O.R.P. #2 

For the second trajectory planning problem, the Legendre-Clebsch condition will be 
met in strengthened form if a > 0, since the Legendre-Clebsch necessary condition [76] is 

iL + - «><> (3.46) 

Defining F(x,y,y’) = f 2 + ey 2 + ay’ 2 , the Weierstrass excess function for this problem 
turns out to be 

E (x,y,y',p) = F(x,y,y') - F(x.y.p) - (y'-p)F p (x,y,p) = a (y'-p) 2 (3.47) 

This is positive if y' is not equal to p. 

The Jacobi's differential equation [76] for this problem turns out to be: 

T1 " = a (E + f y + ffyy)T1 (3.48) 

Here T| is the solution to the linearized Euler-Lagrange equations or the second- variation. If 
the coefficient ^ (e + fy + f fyy) on the right hand side is a slowly varying quantity with 

respect to range, additional analytic results may be obtained using the WKB method [16]. 
Let 

q(T) = ^- (e + fy + f fyy) • (3.49) 

where x is the stretched range x = p. x, p. > 0. Equation (3.48) can be rewritten 

Tf-q(x)Ti=0 (3.50) 


The various derivatives in the above differential equation can now be written as 
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dx dx dx *ax 

d 2 ^ _ 2 d2 n 

dx 2 dx 2 


H 2 


d 2 ^ 

dx 2 


q(x) ti = 0 


with change of variable |i = X* 1 
d 2 !] 


dx 2 


- X q(x) r| = 0 


(3.51) 

(3.52) 

(3.53) 


(3.54) 


The equation (3.54) with a large parameter X is referred to as Liouville's problem [74]. Its 
approximate solution can be obtained using WKB method [16]. 


H - qHJ[C 1 cos (J^ ^/qdx) +C 2 sin (£ ^/q"dx)}, when \ (3.55) 



T] <0/ (Cj exp (J ^/^qd\) + C 2 exp (J .y^’qdx)}, when (q< 


(3.56) 


where Cl, C2 are arbitrary constants. It may be observed from the above solutions that if 
q > 0 in the given interval, there are no conjugate points. Since f > 0 by definition, this 
expression implies that q is always positive if fyy > 0. That is, if the vehicle is traveling in 
a valley or along a saddle, the resulting extremal will not contain any conjugate points. 

In the cases where q(x) changes sign along a given extremal, one has to carry out this 
test numerically. The procedure employed in the previous trajectory planning technique 
may once again be used. It may be verified that the integrand of the performance index 


satisfies the convexity condition given in [69]. As a result, if a given extremal satisfies the 
second order necessary conditions, it affords a global minimum. 

Computational Effort 

3.5.1 Com p utational Effort for O . R.P . # 1 

Studies using a VAX 1 1/750 computer have shown that an extremal requiring about 70 

integration steps consumes between 1. 5-2.1 seconds of CPU time. Given the desired 

initial and final conditions on down-range and cross-range, several Euler solutions need to 

be evaluated to converge on the one satisfying the given boundary conditions. Using the 

computational flow chart in Figure 3.2, between six and seven iterations were found 

adequate to reduce the boundary condition error by an order of magnitude. Clearly, this 

will depend on the given terrain profile. For example, if the terrain gradients f* and f y were 

small such that the boundary condition error depends approximately linearly on the initial 

heading angle error, the number of iterations and boundary condition error may be related 

as using the convergence formula for the method of bisections [62] as. 

In (e yi /ey,) (3 ' 57) 

r In 2 

In the expression (3.57), r is the number of iterations, e^ is the initial boundary condition 
error and ey is the boundary condition error after r iterations. It is important to stress here 

that this relationship does not account for nonlinearities due to terrain profile. Its 
usefulness is limited to generating a first order estimate on the number of iterations required 
to satisfy a specified boundary condition error tolerance. 
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In the present research the total computations for one step integration are about 626 
multiplications and divisions and 240 additions and substructions. 

3.5.2 Computational Effort for O.R.P. #2 

Assuming that the stationary condition (3.37) is integrated using the Euler's method, 
and that a linear interpolation is employed for computing the terrain altitude at various 
down-range and cross-range locations, one can obtain a Formula for the required number of 
operations. In the following, it is assumed that multiplication, division, addition, and 
subtraction each count as one operation. Each integration step is found to require 23 
operations. The method of bisections requires two mathematical operations and one logical 
operation per iteration. Thus, if there are n discretization intervals in the down-range 
direction and r bisections iterations, the present method would require approximately 

m = 23 nr (3.58) 

operations. Assuming that it is desired to decrease the interval of uncertainty by two orders 

of magnitude, one requires about 7 iterations using the method of bisections. Thus, the 

number of operations required for the method of bisections is m = 161n. Note that one can 

* 

reduce the interval of uncertainty by eight orders of magnitude by increasing the number of 
iterations by about four times. 

Next, assuming equal discretization of the down-range and cross-range directions, the 
dynamic programming scheme is found to require 23 mathematical operations per node to 
evaluate the performance index. In a typical computing scheme [12], this will have to be 
evaluated at n 2 - 1 nodes. After these computations, one has to execute 2n -1 logical 
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operations. Thus, the total number of operations required in a full scale discrete dynamic 
programming approach for the present route planning problem is approximately 

m = 23 n 2 + 2n ( 3 - 59 ) 

operations. Thus, to first order, the two methods are equal in terms of computational speed 
whenever n = 7. Whenever the number of discretization intervals exceeds this figure, the 
present approach provides a faster solution. In any case, note that the computational effort 
in the present approach is a linear function of the number of discretization interval, while in 
the dynamic programming approach, the computational effort is quadratic in the number of 
discretization intervals. 

On a rough terrain, one would require a large number of discretization intervals to 
obtain results with sufficient fidelity. In that situation, the advantage of the present 
approach will be even more significant. On the other hand, the dynamic programming 
approach guarantees the optimality of the solution, while in the Euler solution method, this 
has to be verified through second order necessary conditions . 

3.6 Wind Effects 

The magnitude and direction of winds are known to have a substantial effect on the 
performance of most aircraft. This effect is accentuated on a conventional helicopter 
because of the low flight speeds. As mentioned in Section 3.2, the ambient winds can be 
incorporated into the trajectory programming scheme with slight increase in complexity. 
For the purposes of illustration, the effect of constant winds along x-axis in the first 
trajectory planning problem will be considered. Using equations (3.8), (3.10), and (3.1 1), 
and assuming the same performance index (3.14), one can obtain costate equations as 



(3.60) 


- { 1 - K + K f(x,y)} Vl + ft cos x 
V + u Vl + ft cos x 

( 1 - K + K f(x,y))(V 1 + f| + sin x - f x f y cos x) 
Vl + ft (V + u Vl + ft cos x) 


Following the same steps as in Subsection 3.3.1, the differential equation for heading 
angle x turns out to be 

... V (Ci sin x + C 2 cos x) + u Ai cos x Pi sin x + F>2 cos x) 

X x x 

(1 -K + Kf) A 3 A| 


where, 

Ai=Vl+ft 

A 2 = Vl+ ft + fy 

Bj = {- Al f x f xx + Aj (f y f xy + f x f xx )} Ai 

B 2 = - { fxfxx A 2 + Aj(f y f xy + f x fxx) } fxfy + Aj A 2 (f x fxy + fyfxx) 

Ci = - Kf x Aj Aj - (1 - K + Kf) BjA 2 

C 2 = - Kf y A? A 3 2 + (1 - K + Kf) { B j f x f y + B 2 + A&f» (f x f y - f xy A?) } 
Di = - Kf x Aj k\ + (1 - K + Kf) (- BjA 2 + A 3 2 f x f xx ) 

D 2 = - KA]A 3 2 (1 + fy) 2 + (1 - K + Kf) (B |f X fy + B 2 ) 


The accuracy of these expressions may be verified by putting the wind speed u to zero. 
In this case, the equation (3.62) exactly corresponds to (3.26). Figures 3.17 and 3.18 
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show several extremals for the minimum flight criterion and the maximum masking 
criterion when wind to airspeed ratio u/V equals to 0. 1 and all conditions being the same as 
the Figures 3.4 and 3.5. An observation that may be made from these figures is that the 
winds affect the maximum masking trajectory to a higher degree. 

It is possible to include wind effects in the second route planning problem also. 
However, this will not be pursued in the present research. 

3.7 Conclusion 

In this chapter, two systematic methodologies for optimal trajectory planning (ORP) 
useful in the helicopter terrain-following/terrain-avoidance flight was presented. 

For the ORP No.l, the terrain constraint was embedded into state equations via a 
coordinate transformation. The performance index consisted of a linear combination of 
flight time and terrain masking. Using an adjoint-control transformation, the optimal 
control problem solution was reduced to a search for the initial value of heading angle. It 
was shown that the optimal airspeed, a second control variable in the formulation, should 
be chosen as the maximum permissible value. A simple computational scheme based on 
the method of bisections and a fifth order Kutta-Merson numerical integration technique 
was outlined for generating Euler solutions. Families of Euler solutions for minimum 
flight time and maximum terrain masking were presented. 

It was shown that the Legendre-Clebsch necessary condition and Weierstrass excess 
function are satisfied everywhere in the admissible region. Further, conjugate points have 
been shown to occur in certain regions of the specified terrain. In the regions where 
conjugate points do not occur, the Euler solutions for the maximum masking problem 
provide a strong local minimum. The Euler solutions for the minimum time problem 
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satisfying the second-order necessary conditions were shown to provide a global 
minimum. 

For the ORP No.2, the performance index consisted of a quadratic form in the terrain 
altitude, lateral deviation from the nominal trajectory, and heading angle. By changing the 
independent variable from time to down-range, the order of the problem was reduced. 

The Legendre-Clebsch necessary condition and Weierstrass excess function were 
satisfied everywhere in the admissible region. For this problem, an approximate conjugate 
point test was developed using the WKB method. 

The winds effects on the trajectories were briefly examined for ORP No. 1. It is found 
that the winds affect the terrain masking trajectory to a higher degree. 
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CHAPTER IV 

TRAJECTORY PLANNING AS A TWO-SIDED 
OPTIMAL CONTROL PROBLEM 

4.1 Introduction 

Previous chapter dealt with the situation wherein we seek to determine the optimal 
trajectory for a vehicle to fly from an initial position to a given final position. In this 
chapter we will examine the same problem involving two vehicles with conflicting 
objectives. 

In such a situation, the trajectory planning problem will turn out to be a differential 
game. An example is the case where one of the vehicles is attempting to intercept the other 
while flying in a terrain following mode to avoid detection. Meanwhile, the second vehicle 
may be executing a terrain following flight with the objective of avoiding capture. 
Problems of this nature have received scant attention in the previous literature [77]. 

Next generation military helicopters such as the LHX must have capabilities to 
automatically engage with ground and air targets. The avionics requirements for the low 
altitude air combat mission are complex, and reflect the problems brought about by 
simultaneous air combat and terrain flight [78]. 

Since the publication of Isaac’s famous book [17], the homicidal chauffeur problem 
has became a model for vehicle pursuit-evasion. The solutions to this problem are very 
intricate and have been discussed in References 79 and 80. Ciletti [77] indicated that the 
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assumption of perfect information and role ambiguity [81] are some of the factors that have 
prevented the application of differential games theory in realistic military situations. 
Additionally, dynamics of the vehicles are highly nonlinear and of high order. These 
factors have led to a situation wherein complete results are given only for a very limited 
class of problems. 

In this chapter, two different approaches are proposed for studying the helicopter 
pursuit-evasion problem. A nonlinear pursuit-evasion game employing helicopter model 
incorporating the terrain profile is discussed first. This model was used in Chapter in for 
trajectory planning. This is followed by an approach based on the feedback linearization 
technique motivated by the research in Reference 91. Each of these approaches will be 
discussed in the following sections. In the ensuing, it will be assumed that each player has 
complete information on the helicopter parameters and a noise free measurement of all the 
state variables. 


4.2 A Nonlinear Pursuit-Evasion Game 
4.2.1 Probl em Formulation 

The helicopter pursuit-evasion problem is analyzed here as a deterministic two-person 
zero-sum differential game. The game begins at a certain set of initial conditions, when the 
helicopters first become aware of each other. In the present analysis the respective roles of 
the players are assigned at the outset. It is assumed that this role definition remains 
unchanged during the entire game. Although this formalism restricts the applicability of the 
results, it is useful in revealing salient solution features. In real situations, the evader may 
be a helicopter without offensive capabilities or with limited maneuvering capabilities. 
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4.2.1. 1 Equ ations of Motion 

When compared with high performance aircraft pursuit-evasion, one of the 
characteristics of helicopter pursuit-evasion [82, 83] is that each participant attempts to 
reduce its exposure to the opponent by avoiding high altitude flight throughout the 
engagement. This is due to the fact that the vertical maneuvering above NOE altitudes may 
make the helicopters more vulnerable to detection by ground and air based surveillance 
system. The equations of motion developed in Chapter in are used here to model the 
helicopter flight in terrain-following mode. In the interest of clarity, these equations of 
motion are repeated here. 


Ve COS %e | Ye fx.fy.sin % e 

Vi+f£ V i+f £V i+f£+fy. 

^ = YeVl-t-fj sin %e 
y ' = ' 

t Vp cos Xp t Vpfx p fy p sinXp 
P VT< V l+fi p V l+fx P +fy P 

VpVl+f2 p sin Xp 


(4.1) 


(4.2) 


(4.3) 


(4.4) 


The subscript e denotes the evader, while the subscript p denotes the pursuer. To make the 
problem meaningful, it is assumed that the pursuer has a higher speed than the evader, i.e., 


Vp>V e 


(4.5) 
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Whenever the terrain gradients f x , f y are zero, it is possible to show that the condition (4.5) 
is essential to guarantee capture. The initial conditions at t = 0 are 


(x«,ye)lt=o specified 

(4.6) 

(x p ,y p )l t =o specified 

(4.7) 


This problem is a terrain following version of simple motion illustrated graphically in 
Reference 17. The participants select the direction of travel and may change it as fast as 
desired. Thus, the control variables in this problem are the heading angles, % p and Xe- 
Note that this pursuit-evasion game has a very simple solution whenever the terrain 
gradients are zero. However, the solution becomes quite complex in the general case; as 
will be apparent in the ensuing. 


4.2.1.2 Termination of Game 


The game terminates at the first instant the pursuer succeeds in approaching the evader 
within the firing range of its weapon system. For the case of a circular weapon range 
envelope, this condition can be expressed by the requiring that 


(x e - x p ) 2 + (y e - y p ) 2 Itt^d 2 
4-{(Xe-x p ) 2 + (y c -y p ) 2 } <0 


(4.8a) 

(4.8b) 


The quantity d is a specified number. The condition (4.8b) is required to ensure sufficient 
time for weapon usage. The time of capture is determined from the requirement that 

t f =min {(x e -x p ) 2 +(y e -yp) 2 -d 2 =0} 


(4.9) 


r 


65 


i.e., the instance that the condition (4.8a) is met as an equality. 

It has been suggested in the literature [17] [84] that a fan-shaped capture set may be 
more realistic. However, in all that follows, a circular envelope will only be considered. 

4.2. 1.3 Measure of E ffectiveness 

In the most planar games reported in the literature [80] [85], the performance index of 
the game is the time of capture. The objective of the pursuer is to minimize the terminal 
time tf, while the evader endeavors to maximize it. Additionally, in the present setting both 

players also attempt to minimize their flight altitude to ensure adequate terrain masking. In 
this case, the performance index can be expressed as 

‘r 

J = min max J(l+ W p f p - W.fjdt (4.10) 

*F * • t 

0 

where, W p and W e are weighting factors for the pursuer's altitude and evader's altitude, 
respectively. The negative sign on the second term explicitly recognizes the fact that the 
evader is attempting to maximize the performance index. In order to satisfy the terminal 
constraint (4.8), it may be appended to the performance index in the form of a penalty 
function [12]. In this case, the augmented performance index is 

J a = min max Q(t f ) + J (1 + W f p — W e fe) dt (4.11) 

X P x. t. 


The terminal penalty function is given by 

Q(tf) = ^ { (x p - x*) 2 + (y p - y c ) 2 - d 2 } l t( (4. 1 2) 


and v is an undetermined multiplier in the game. 
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4,2.2 Derivation of Optimal Strategies 
For convenience, define 


ri= * e . r 2 s ye. r 3 s i p , r 4 = y p 


(4.13) 


The variational Hamiltonian [12] may be formulated as 

4 

H=H-Wpfp-W e f c + X^i 

i=l 

Expression (4.14) may be written in a more succinct form as: 
H = 1 + H(x p ) + H(Xe) 


(4.14) 


(4.15) 


with 


H(Xp) = W p f p + X 3 r 3 + X 4 r 4 


(4.16) 


H(Xe) = -W c fe + Xiri+X 2 r 2 


(4.17) 


Since H is separable in terms of pursuer-evader controls Xp and x e , the saddle point 

condition [39] is satisfied. As a result, the order of maximization and minimization does 
not influence the outcome of the game. Thus, one has 


min max H = max min H 

I, I. I. z F 


(4.18) 


The integrand of performance index (4.10) and system equations (4.1 - 4.4) are continuous 
and satisfy the Lipschitz condition. Thus, the sufficient condition for the game to have the 
value is also satisfied, see Appendix C for more details. Each player's goal is to reach the 
value, i.e., saddle point of the performance index. 
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Euler-Lagrange equations for the evader are given by 

Bi sin Xe ^2 + (B 2 sin Xe + B 3 cos Xe) h 


Wef^- 


A]A\ 


V e 


>-2= W c f ye 


B 4 sin Xe ^2 + (B 5 Sin Xe + B 6 cos Xe) h. y 


a^a 3 2 


(4.19) 

(4.20) 


where 

Ai-Vl+^ 

A2 = Vl+ W e + ly e 

Bi = {- A 2 fx e fxex e + Aj (fy c fxey e + fx c fx e x e ))Ai 

B 2 = - { ^Xe^XeXe A 2 + Ai(fy e f Xe y e + fxe^XeXe) }fx e fy e 

+ A^ A 2 (fx e fxey e + fye^*e*e) 

^3 = " A 2 fx e fx c x e 

B4 = {* A 2 fxe^Xeye + A? (fx c fx e y e + A i 

B 5 = - {fxe^Xey c A 2 + Ai(f Xe f Xe y e + Wyeye^ksfye 
+ A* A 2 (fy e f Xe y e + fxefy c y e ) 

®6 = ' A 2 fx e fxcy e 

The Euler-Lagrange equations for the pursuer can be obtained as 

; „ T , Di sin x P X 4 + (D 2 sin x P + D 3 cos x P ) h „ 

*3 = - W p f Xp 2 — V p (4.2!) 
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_ v,r r D 4 sinXpX4 + (D 5 sinXp + D 6 cosXp)^3 „ 
m - - w pty P * v p 


(4.22) 


Ci = V 1 + f? p 
C 2 = V 1 + f| p + f$ p 

Di = {- c| fxpfxpXp + C? (fy p fx p y p + fxpfx p xp)}c? 

^2 = - {fxpfxpXpC^ + Cl (fy p fx p yp + fxpfxpXp) } ^Xpfyp 
+ Cl (fxpfxpyp + (yp^XpXp) 
f^3 = “ C| ^pfxpXp 

D 4 = (• C2 fxpfx p y p + Ci (fxpfx p y p + fy p fy p y p ) } C? 

^5 = ■ (fxp^XpypCl + Cl (f X pfxpy p + fypfy p y p ) ) fx p fy p 
+ Cl C2 (fy p fx p y p + fx p fy p y p ) 

= * C 2 fxpfxpyp 


The optimality conditions for the two participants are given by 


tan Xe = 


Xlfx e fy e - hj (1 + fj) 

^1 V 1 + fxe + fy c 


tan X P = 


^fxpfyp ~ ^.4 (1 + f| p ) 

*3 V 1 + f| p + f$ p 


(4.23) 


(4.24) 


Terminal conditions on costates can be obtained from the terminal penalty function as 
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Xf(tf) =^-=' v (X p -X c ) 

(4.25) 

fcaOfl = I7- = - v (y P -ye) 

dye 

(4.26) 

M tf )=|^= V(Xp-Xe) 

(4.27) 

^ (tf ) = !r = v (yp-y*) 

dy p 

(4.28) 

a 

II 

sis 

II 

O 

(4.29) 


Since the final time tf is unspecified and the Hamiltonian does not explicitly depend on t, 
this problem has a constant of motion, viz., 

H(t) = 0 , 0 < t <£ t f (4.30) 


Substituting equations (4.25) - (4.28) into (4.14) and invoking (4.29), it is possible to 
obtain an expression for the undetermined multiplier v as: 

v- 1+Wpf p -W e fe I (4.31) 

(Xp-Xe)ri+(yp-ye)r2-(xp-Xe)r 3 -(yp-ye)r4 f 


It may be verified that the denominator of equation (4.31) is simply the negative of the 
product of range and range rate. Since the final value of the range is positive and the 
terminal range rate is negative, the parameter v has a positive sign for the pure pursuit- 
evasion game in which Wp = We = 0. Moreover, expression (4.31) implies that at the final 
time, if one requires the pursuer and evader’s positions to match exactly, v would become 


infinite. 
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4.2.3 Numerical Results 

Most differential game solutions reported in the literature are obtained using numerical 
methods. Among them, simple shooting in retrogressive time [86] [87], differential 
dynamic programming [88], and gradient method [89] have been widely used. In certain 
situations, the solutions obtained using reduced order modeling can be corrected for 
neglected dynamics using singular perturbation techniques [90]. 

In the present research, retrogressive time integration is used to generate the extremals. 
The trajectories are generated by first selecting the terminal position of the evader. The 
terminal position of the pursuer is then selected from the specified capture set. This 
corresponds to the expression (4.8a) being met as an equality. Substituting the final 
costates into the optimality conditions (4.23) and (4.24) and noting that the unspecified 
multiplier v exists in both numerator and denominator of these equations, the final values of 
pursuer and evader's heading angle can be calculated. The final values of ri, r 2 , r 3 , and T 4 
can be computed. Substituting these in equation (4.31) yields the unspecified multiplier v. 
This enables the evaluation of the final value of the costates. The state-costate system is 
then integrated backwards in time to obtain the pusuer-evader trajectories. 

Figure 4.1 shows the pursuer and evader trajectories with terrain masking weights set 
to zero. As in the one-sided trajectory optimization problem, the trajectories are nearly 
straight lines. In Figure 4.2, the pursuit-evasion trajectories with Wp = 1, We = 1 are 
illustrated. The value of the parameter v is computed as v = 5.79E-3. With nonzero Wp 
and We, the trajectories exhibit significant terrain masking features. Figure 4.3 illustrates 
the altitude evolution as a function of time-to-go. Note that tg 0 = 0 corresponds to the game 
termination instant. Pursuit-evasion trajectories for another set of initial and final 
conditions with Wp = 1, We = 1 are given in Figure 4.4. The value of the parameter v is 



computed in this case as v = 2.03E-3. The corresponding altitude history is given in 
Figure 4.5. Figure 4.6 illustrates trajectories for the pursuer and evader with Wp = 0.0, 
We = 1.0. Note that in this case, the pursuer does not attempt any terrain masking. As a 
result, its trajectory is nearly a straight line joining the initial and final conditions. The 
value of the parameter v is computed as v = -.747E-3. For Wp = 0.5 and We =1.0, Figure 
4.7 shows the pursuit-evasion trajectories with v = -.375E-3. For the same weighting We, 
the trajectories become more curved as the weighting Wp is increased. 

Legendre-Clebsch Condition 

The optimality conditions (4.23-4.24) do not guarantee that the players maximize or 
minimize the performance index according to their role. To determine the correct optimal 
controls that make the pursuer minimize the performance index, while making the evader 
maximize, the verification of Legendre-Clebsch condition is necessary . 

For the present differential game, the Legendre-Clebsch condition requires that 

^ 0 (4.32) 

HjmCp ^ 0 (4.33) 


where 


Hx*. “ 




COS Xe+ (^-2 (1 +*2 e )- klWyJ sin y 

Vl +fl* V 1 + fx« + fy e 

- X 3 V 1 + f% p + f$ p cos x P + (Ml + ftp) - Vxpfyp} sin % p ^ 

Vi + f| p V 1 + f| p + f$ p 


(4.34) 

(4.35) 


The conditions (4.32),(4.33) may be strengthened by requiring strict inequality. 
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Since these expressions are too difficult for hand computations, they are numerically 
evaluated along the extremals. Figure 4.8 shows these quantities evaluated along the 
trajectories in Figure 4.2. It may be seen that the Legendre-Clebsch necessary conditions 
are satisfied in the strengthened form. Figure 4.9 illustrates the L-C test along the 
extremals in Figure 4.4. In this case, it appears that this test is satisfied in the strengthened 
form every where except a small region near the final region. It is important to emphasize 
that the Jacobi test needs to be carried out along these trajectories before concluding the 
optimality. This will not be pursued any further in this thesis. 


4.3 Feedback Linearized Solution to a Pursuit-Evasion Game 

As illustrated in the previous section, the numerical solution of even simple differential 
games requires a tremendous amount of computational effort An approach proposed 
recently for solving a class of differential games [18] is to transform a nonlinear model into 
linear time-invariant form. A linear differential game is then solved and the results are 
transformed back to original coordinates to obtain a nonlinear feedback law suitable for real 
time implementation. The advantage of this approach is that it can handle high-order 
nonlinear vehicle models In the analysis. In the following, such a differential game will be 
formulated and solved. The ensuing formulation will employ a high order model of the 
helicopter. 

4.3.1 Vehi cle Model 


The point-mass model for a high performance helicopter can be expressed by the following 
nonlinear differential equations: 
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= Vcosycos y 

(4.36) 

y = V cos y sin y 

(4.37) 

h = V sin y 

(4.38) 

V = T sin 0 . g sin y 

(4.39) 

. g Tcos0cos<l> 

Y=V< mg •'**’/> 

(4.40) 

• T cos 0 sin 6 

\1 f = 

mV cos Y 

(4.41) 


Here, x is the down-range, y the cross-range, h altitude, V the speed, y the flight path 
angle, y the heading angle, T the main rotor thrust, m the vehicle mass and g is the 
acceleration due to gravity. Figure 4.10 shows the definition of the axis system. The 
control variables are the pitch attitude of the helicopter 0, roll attitude of the helicopter <>, 
and the main rotor thrust T. If desired, one may model T in terms of two components. 
The first depending only on the vehicle states To and the second component 5T that 
depends purely on the collective control, i.e., T = To + 8T. The assumptions involved in 
this model are stationary atmosphere and uniform gravitational acceleration. It is important 
to note that this highly simplified helicopter model does not permit hovering. 

4.3.2 Problem Formulation 

The point mass helicopter model (4.36) - (4.41) will be transformed to a linear, time- 
invariant form in this section. Various steps involved in this transformation will be 
discussed in detail. 


In order to execute the terrain flight mode, admissible trajectories should satisfy the 
terrain profile constraint i.e., the vehicle altitude h should be: 


h = he + f(x,y) 


(4.42) 


This equality constraint can be absorbed into the state equations as will be shown in the 
following. Differentiating the equation (4.42) with respect to time and noting that h c is a 
constant, one has 

h = fx* + f y ? ( 4 . 43 ) 

Equating (4.43) to (4.38), and substituting (4.36) and (4.37) for X and $ results in 

tany=f x cosy + f y siny (4.44) 

Differentiating once again with respect to time results in, 

scc^yy = f xx X cosy + f xy y cosy+ f^x siny + fyyy siny - f x siny y + f y cosy y 

(4.45) 


Next, substituting from (4.36), (4.37), (4.40), and (4.41) for X, y, yandy, and 
rearranging yields the following relation 

* 

T cos0 cos<t> = mV^cos^ { (f w cosy + f xy siny) cosy + (f xy cosy + fyy siny) siny } 

+ (- f x siny + f y cosy) cosy T cos0 sin<j> + mg cosy (4 4 ^ 

The expression (4.46) relates the vehicle states and control variables to the terrain profile. 

It specifies the vertical force required to ensure that the vehicle trajectory satisfies the 


terrain profile constraint. As a result, this expression dictates the altitude dynamics for the 
two vehicles. 
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Following Reference 18, the remaining components of the nonlinear helicopter models 
for the pursuer and evader may be transformed into the Brunovsky canonical form by 
differentiating the expressions (4.36) and (4.37) with respect to time and substituting for 
V, y, \j / from the expressions (4.39) - (4.41). Defining four pseudo control variables, the 

helicopter models can be put in the form 


xp ” <*pl 

>T3 

II 

& 

(4.47) 

II 
1— » 

y e — ac2 

(4.48) 


where, a p i, 5 p 2 , aei and ^2 are the vehicle acceleration components in the earth-fixed 
frame. These quantities may be related to the pursuer-evader states and controls as 


T„ sin 0 


Xrj am Un . v 

a p i = ( p m p - g Sin Yp) COS Y P cos \|f p 


T p cos 0 P sin 4> p sin y p 


m r 


T„ cos0 D cosd>p . . 

g (-2 — g^g — — - cos Y P ) sin Y P cos \|r p 


(4.49) 


.Tp Sin 0n . • 

a p2 = ( p m g - g sin y p ) cos Y P sm y p 


T p cos 0 P sin 4>p cos y p 
+ m„ 


T D cos 0 D cos <b D . . 

’ S (-* — ^ - cos Yp) sm Yp sin Vp 


(4.50) 


T sin 0 T e cos 0 e sin <J> e sin 

ael = (iijgiSt. - g sin Ye) COS Ye cos Ve - 


-g(- 


me 

Te COS 0 e COS <t> e 


me 


nieg 


- COS Ye) sin Ye COS \|/e 


(4.51) 


T sin 0 v . T e cos 0 e sin <|>e cos y e 

5a = - g sin Y«) cos Y, Sin — 


-g(- 


me 

Te COS 0e cos <t> e 


meg 


- cos Ye) sin Ye sin \jr e 


(4.52) 
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Equations (4.47) and (4.48) are linear time-invariant form. If the pseudo control variables 
5 p i, S p 2 , Sei and ^2 were known, the actual helicopter controls can be obtained using the 
following equations: 


< D p = tan-i(V COS¥ P' I P 1 $in yP.) 


0 = , sin <|>p (api cos y p + a p2 sin y p + F p sin y P ) 

p (Sp 2 cos v P - 2pi sin y p ) cos y p 


2 p 2 cos \|/p - ap! sin Y p 

l p ~ a . , 

cos 0 P sm 4> p 


“p 


(4.53) 

(4.54) 

(4.55) 


where, F p = (a p 2 cos Y P - a p i sin Y p ) cos y P (-f* p sin Yp + fy p cos Y P ) + g cos y p 
+ Vj cos 3 y p (f XXp cos 2 Y P + 2 f yXp sin Ypcos Y P + fyy p sin 2 y p ) 


Note that the equations (4.53) - (4.55) used the vertical force component constraint 
expressed in (4.46). A similar set of expressions may be obtained for the evader also. 


4.3.3 Linear Quadratic Pursuit-Evasion Game 
With the definition of position vector 

(4.56) 

the feedback linearized vehicle models (4.47), (4.48) may be expressed in the standard 
state variable form as: 

f p*V r pfto) given (4.57a) 



J 


t e = v e . r e (to) given 


(4.57b) 



1 
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Vp = £tp , Vp(to) given 

(4.57c) 

- 

v c = a e , v e (to) given 

(4.57d) 

- 

Next, introducing the relative coordinates, one has 



r = r P -r e 

(4.58a) 


< 

It 

i 

< 

o 

(4.58b) 


These dynamic equations can be written in a more succinct form as 

CHJMII-Ch-* 

(4.59) 


Where 0 is a 2 x 2 zero matrix and I is a 2 x 2 identity matrix. The control variables in the 
model are the acceleration components in the earth fixed coordinate system. The pursuer 
uses the control a p to attempt to capture the evader, while the evader uses the control a« to 
avoid capture. 

To further simplify notation, write the above equation as 


% = F £ + G u 

(4.60) 


with 



>-[ ; j] «-(!) w:i -v- 

' 

- 

The objective of the pursuer is to minimize the terminal miss, which the pursuer 
attempts to maximize. The terminal miss is defined here as 


2 rT (*f)S{T(tf) 

(4.61) 
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The superscript T denotes the transpose of the vector and Sfis a positive semidefinite 
matrix. Integral quadratic acceleration constraints are next imposed on both the pursuer and 
the evader to make this game meaningful, and also to enable the application of the well 
known Linear-Quadratic game results [12]. These constraints are included using two 
positive definite weighting matrices W p and W e . 

Adjoining these constraints to the performance index, one has 

1 1 tf 

J = min max y r T (t f )S f r(t f )+y J (ajw p a p -ajw e a e )dt (4 52 ) 

a p a « t 0 


The saddle-point solution to this problem is given by [12] 

^ = -W p ‘G T s[J] 

a. = -W;'G T s['] 


The matrix S is the solution of the matrix Riccati equation 
$ = -SF - F t S + SGfWp' 1 - W e '*)G T S 


(4.65) 


with the terminal condition 


S(t f ) = 


S f O' 

0 0 . 


In order to obtain solutions in the general case, the matrix Riccati equation (4.65) has 
to be integrated backwards in time. To illustrate the present solution further, 
simplifications will be introduced in this problem. 

If Sf, W p , and W e were constant diagonal matrices with Sf(i,i) = cr*, W p (i,i) = 1/ctj 
and W e (i,i) = 1/Pi, i = 1, 2, the saddle point solution can be expressed in the form 
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ap = kj r + k 2 v 

(4.66) 

a« = k 3 a p 

(4.67) 

In the expressions (4.66) and (4.67), ki, k 2 and k 3 are 2 x 

2 diagonal feedback gain 

matrices with 


k,(M) = — ,1-1,2 


-J- + («i - PiHf 

Oj J 

(4.68) 

k 2 (i,i) = ki(i,i) t g0 , i = 1, 2 

(4.69) 

k 3 (U)= I 1 , i- 1.2 

(4.70) 


The time-try go required for the computation of feedback gains is defined as 

tgo = tf - 1 (4.7 1 ) 

While implementing this guidance law, it is preferable to compute t g0 using feedback 
information. This is because the vehicle model is approximate. Consequently, the 
specified final time will not be equal to the actual time for capture. In the following, a 
method for computing the time-to-go will be discussed. 

4.3.4 Term inal Time Estimation 

Guidance scheme discussed in the foregoing requires an accurate estimate of time-to- 
go for satisfactory operation. Several methods for calculating this quantity are available in 
the missile guidance literature [91, 92]. However, most of the reported time-to-go 
calculation methods neglect the two-player nature of an engagement scenario. In Reference 
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93, an exact method for calculating time-to-go was outlined. This approach is developed 
for the present differential game in the ensuing. 

Since the transformed min-max problem in Subsection 4.3.2 does not explicitly 
depend on time, the variational Hamiltonian is constant, i.e., 

H(to) = H (t) = H(tf) = constant (4.72) 

If the desired conditions at the termination of the differential game can be defined in terms 
of the final differential position and velocity components, it is possible to evaluate the 
variational Hamiltonian at the final time. This is feasible in the present problem since the 
control variables are available in state feedback form with time-to-go as the only unknown 
parameter. Next, equating the numerical value of the corresponding expression at the initial 
time results in a polynomial in time-to-go. A positive real solution of this polynomial is 
then the exact value of time-to-go. In case of multiple solutions, the smallest value may be 
selected. For further details on the calculation of time-to-go, refer to Reference 93. 

For the present application, assuming that the weighting matrices Sf, W p , and W e are 
identity matrices multiplied by scalars a, 1/a, and 1/p, respectively, the polynomial 
equation for the terminal time is as follows [93]: 

R 2 { £ + — ) 2 = (Ax + Ax t f ) 2 + (Ay + Ay tf) 2 ( 4 . 73 ) 


where 


R 2 = {(x p - x«) 2 + (y p - y e ) 2 }ltt 


Ax = (x p - x e )lto , A* = (x p - 


Ay = (y P - y e )lto . Ay = (y p - y e )ito 



Note that the quantity R defines the capture radius. 

If the current time is assumed to be the initial time and set to zero, t g0 = tf. The roots 
of this sixth order polynomial are the values of time-to-go. Since there are six possible 
values, the smallest real positive value needs to be used. In the present work, the roots 
were found using the method of golden section. 

4.3.5 Numerical Results 

Two scenarios are examined for the pursuer and evader starting at the coordinates 
(4000, 4000) and (5000, 5000), respectively. In both engagements, the weighting factors 
are chosen as follows: cr = 1.0, a = 0.0007, and (3 = 0.00001. 

Figure 4.11 shows the trajectories for the pursuer and evader both with zero initial 
heading angles. The corresponding velocity histories are shown in Figure 4.12. The 
speed variations occurring due to the terrain profile can be seen to affect the pursuer to a 
higher degree. This is because the pursuer is assumed to have a higher acceleration 
capability when compared with the evader. Load factor histories for two helicopters are 
given in Figure 4.13. The load factor is defined here as the ratio of main rotor thrust and 
the helicopter weight. Evader's thrust appears to be smoother than the pursuer's. Figures 
4.14 and 4.15 illustrate roll and pitch attitudes of helicopters. Various features appearing in 
these figures arise from the terrain profiles. Altitude histories for two helicopters are 
shown in Figure 4.16. 

In the the second scenario the pursuer has an initial heading angle 90 angles and the 
evader is at a 0 degree heading angle. The resulting trajectories are shown in Figure 4.17. 
As may be observed in Figure 4.18, this game involves the typical tum-and-dash strategy. 
To tum quickly, the pursuer first decelerates and then accelerates to catch the evader. 
Figure 4.19 shows load factor histories. Figures 4.20 and 4.21 show the corresponding 
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roll and pitch attitude histories. Altitude histories for two helicopters are given in Figure 
4 . 22 . 


Conclusion 

The helicopter pursuit-evasion problem while executing the terrain-following/terrain- 
avoidance mode was studied in this chapter. Numerical methods for obtaining solutions to 
these problems were outlined. As an alternative to numerical method, feedback linearizing 
transformations were combined with the linear quadratic game results to synthesize explicit 
nonlinear feedback strategies for helicopter pursuit-evasion. Further investigation of these 
solutions will be of future interest. 



CHAPTER V 


PERFORMANCE VERIFICATION 


5.1 Introduction 

The one-sided optimal trajectory planning schemes of the helicopter discussed in 
Chapter in considered only the kinematic equations. The objective there was to make the 
optimal control problem analytically tractable. To verify whether the generated trajectories 
satisfy the physical constraints, these need to be evaluated on a detailed simulation of a 
helicopter . 

In this report, a six degree of freedom helicopter simulation program called the 
"TMAN" [21] is used for evaluating the generated trajectory. This program was developed 
from a more general helicopter simulation program called ARM COP [94]. This computer 
code is being used in the Vertical Motion Simulator to conduct piloted simulation of 
helicopter flight including air combat at NOE flight levels. Previously, a path planning 
scheme using dynamic programming [8] has been tested at the NASA Ames Research 
Center using this program. 

The TMAN program incorporates a six-degree-of-freedom rigid body vehicle model, a 
first order engine lag response, a simple closed form trim solution, and linearized quasi- 
static aerodynamic force and moment equations. Originally, this program was designed to 



84 


simulate a generic helicopter motion in response to a joy stick. As a result, a coarse 
stability augmentation system is incorporated in the simulation. Recently, Heiges [19, 59] 
developed a trajectory controller for this helicopter model using feedback linearization. The 
TMAN program together with the controller developed in [59] is used in this report to 
evaluate some of the trajectories generated in Chapter HI. 


5.2 Simulation Results 

Trajectories given in Figure 3.6 are used for evaluating the helicopter performance. 
Note that these trajectories were generated with the constant airspeed assumption. The 
airspeed employed in the present investigation is extracted from the available helicopter 
performance data. For example, Reference 82 has presented AH- IS data for both 
maximum rate of climb and rate of descent over its entire speed regime. This data is 
summarized in Table 5.1. It may be observed from this table that over the density altitude 
1000 to 3000 feet range, the maximum rate of climb varies between 8 ft/sec to 36 ft/sec. 
The best rate of climb appears to occur at an airspeed about 100 ft/sec. Since the TMAN 
program uses an AH- IS type helicopter data, all simulations were carried out at this speed. 

Figures 5.1 through 5.8 illustrate the helicopter responses for the maximum terrain 
masking trajectory while Figures 5.9 - 5.16 give the simulation results along the minimum 
flight time trajectory. Figures 5.1 and 5.9 show the altitude rate for both trajectories. The 
maximum altitude rates are about 35 ft/sec. The nonlinear controller [19] was designed on 
an assumption that the cyclic stick and pedals are primarily moment generating controls and 
do not make a significant contribution to the body forces. The collective is the only direct 
force control. Since airspeed is constant in these simulations, altitude change primarily 
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affects the collective control. Figures 5.5 and 5.13 give collective control histories for both 
trajectories. 

As pointed out in Reference 2, the flight along the maximum terrain masking requires 
frequent and severe rolling and yawing motion than the flight along the minimum time 
trajectory. This is because the helicopter seeks to fly at low altitudes to maximize terrain 
masking. The maximum bank angles for the terrain masking and minimum time trajectories 
are about 12 and 6 degrees as shown in Figures 5.7 and 5.15, respectively. The maximum 
heading angle change in the case of maximum masking trajectory is nearly 40 degrees as 
depicted in Figure 5.8. The maximum heading change in the case of minimum time is only 
15 degrees (see Figure 5.16). From Figures 5.6 and 5.14, it is evident that the minimum 
pitch angle is approximately -6.0 for both cases. 

These simulations indicate that the synthesized trajectories are implementable in 
helicopter simulations provided that the airspeeds along these paths are chosen carefully. It 
may be necessary to synthesize an acceleration control loop in the helicopter model to 
ensure adequate ride quality. Clearly, the man-machine interface issues have to be 
sufficiently addressed before attempting a full scale piloted simulation. These issues will 
be of future research interest. 



CHAPTER VI 


CONCLUSIONS AND FUTURE RESEARCH 


This report systematically analyzed trajectory planning schemes for the helicopter 
terrain-folio wing/terrain - avoidance flight by employing optimal control theory and 
differential game theory. Numerical algorithms for trajectory synthesis have been 
developed and validated through simulation. With adequate computing resources, 
trajectory planning methods developed here appear to be implementable on-board the 
helicopter. 


6.1 Conclu ding Remarks 

For the first optimal route planning method, the terrain constraint was embedded into 
state equations via a coordinate transformation. The performance index here consisted of a 
linear combination of flight time and terrain masking.- Using an adjoint-control 
transformation, the optimal control problem solution was reduced to a search for the initial 
value of heading angle. It was shown that the optimal airspeed, a second control variable 
in the formulation, should be chosen as the maximum permissible value. A simple 
computational scheme based on the method of bisections and a fifth order Kutta-Merson 
numerical integration technique was outlined for generating Euler solutions. Families of 
Euler solutions for minimum flight time and maximum terrain masking were presented. 
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It was shown that the Legendre- Clebsch necessary condition and Weierstrass excess 
function are satisfied everywhere in the admissible region. Further, conjugate points have 
been shown to occur in certain regions of the specified terrain. In the regions where 
conjugate points do not occur, the Euler solutions for the maximum masking problem 
provide a strong local minimum. The Euler solutions for the minimum time problem 
satisfying the second-order necessary conditions were shown to provide a global 
minimum. 

For the second optimal route planning scheme, the performance index consisted of a 
quadratic form in the terrain altitude, lateral deviation from a nominal trajectory, and the 
heading angle. By changing the independent variable from time to down-range, the order 
of the problem was reduced. Once again, die Euler solutions for this problem was shown 
to require a one dimensional search. The numerical flow chart developed for the previous 
trajectory planning scheme was shown to be useful for this problem. 

The Legendre-Clebsch necessary condition and Weierstrass excess function were 
shown to be satisfied everywhere in the admissible region. For this problem, an 
approximate conjugate point test was developed using the WKB method 

The wind effects on the trajectories were briefly examined for the first optimal route 
planning problem. It was found that the winds affect the terrain masking trajectory to a 
higher degree than the time optimal path. 

Trajectory planning problem was next formulated as a differential game to synthesize 
optimal trajectories in the presence of an actively maneuvering adversary. Numerical 
methods for obtaining solutions to these problems were outlined. As an alternative to 
numerical method, feedback linearizing transformations were combined with the linear 
quadratic game results to synthesize explicit nonlinear feedback strategies for helicopter 
pursuit-evasion. 
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The one-sided trajectory planning schemes were based on the kinematic equations of 
the helicopter. To verify whether the generated trajectories satisfy the physical constraints, 
these trajectories were tested on a 6 DOF helicopter simulation using a currently available 
flight path controller. The results indicate that the synthesized trajectories are 
implementable provided the airspeed along these paths is chosen carefully. 


Future Research 

Based on the results of this study, the following research areas are recommended for 
further research into the NOE guidance problem. 

Real Time Simulati on and Flight Testing 

As mentioned in Chapter V, NASA Ames Research Center is testing a path planning 
scheme based on the dynamic programming method. The present trajectory planning 
scheme appears to be a viable alternative candidate, thus it needs to be evaluated on the 
Vertical Motion Simulator. Such an investigation would reveal various issues relating to 
both man-machine interface and helicopter physical constraints. 

Use of Alternative Performa nce Indices 

Use of performance indices other than minimum time and maximum masking need to 
be investigated. For example, a weighted combination of time, masking, and flight path 
angle can be considered. Constraining the flight-path angle would have the effect of 
limiting the helicopter climb/descent rate. Since the rate of climb for a helicopter is related 
to the excess power, such a limit on the climb rate may be required in real applications. 
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Helicopter Dynamics in NOE Flight 

The NOE operation requires helicopter fly close to the ground to avoid detection and 
operate at low speeds to avoid collisions with unknown obstacles. This low-altitude/low- 
airspeed flight regime is unique to the helicopter, and yet the least studied out of all the 
flight regimes. Operation in the close proximity of ground brings in the need for studying 
this flight regime in greater detail. 

Time-varvi ng Obstacl e Avoidance 

Known and stationary obstacles may be incorporated into the trajectory planning 
problem by overlaying the artificial envelope on the terrain map as discussed in Chapter HI. 
However, there exist only few path planning schemes examining the problem of avoiding 
the time-varying obstacles, such as schemes for thunder-storm avoidance. References 95 
and 96 may be useful in understanding how the time-varying obstacle can be included in 
the analysis. 


These and other research items will be of future interest. 
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APPENDIX A 


TANGENT PLANE AND THE COORDINATE TRANSFORMATION 


A.1 Introduction 

The purpose of this appendix is to define the tangent plane oriented coordinate system 
used in the development of the route planning problem discussed in Chapters HI and IV. 
Further, the development of a transformation relating vectorial quantities in the tangent 
plane coordinate system to an inertial frame will be developed. 


A-2 Tangent Plane 

Referring to Figure 3.1, let P(xo, y 0 , Zq) be any point on the surface z = f(x,y). If 
f(x,y) is differentiable at (x<>, y 0 ) then the surface has a tangent plane at P. If f x and f y are 
the gradients of the surface along x and y directions, then the equation describing this plane 
is given by [74] 

-f X (x 0 , yo) (x - Xq) -fy(x 0 , yo) (y - yo) + (z - z<x) = o (A. 1) 

Moreover, the vector 

N = - f x l - f y T +k (A.2) 

is normal to the surface at point P. A line which is normal to the surface at P has 
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parametric equations 

x = Xo - f x ((x 0 , y Q ) t 

y = y 0 - f y ((x 0 , y 0 ) t (A.3) 

Z=Zq+ t 

where t is an arbitrary parameter. 


A.3 Coordinate Tr ansformation 

The following discusses the transformation of quantities in tangent plane coordinate system 
to the inertial frame. From (A.2), unit normal vector is given by 


n = 


~ f x i ~ fy j 

V 1 + f x 2 + fy 2 


(A.4) 


Select a unit tangent vector as 


ti = 


i±M 

VTTf 2 


(A.5) 


This unit tangent vector direction is parallel to the plane containing the inertial coordinates x 
and z. Next, the cross-product between the tangent and normal vectors may be used to 
obtain a unit vector orthogonal to these two vectors using the cross-product as follows: 


t 2 = ti x n 


(A.6) 


i.e., 

f X fyl-(l+fg) j-f y k 

Vi + fxV i + $ + 


(A.7) 
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On the tangent plane, the kinematic equations of motion of a particle can be described 


as 


*t = V cos x 
y ( = V sin % 


(A. 8) 


Here, % is the angle made by the velocity vector V with the tangent vector ti. Note that the 
component of the velocity vector normal to the surface is zero. This is because of the fact 
that the vehicle is executing a terrain following flight. The subscript £ denotes the local 
coordinate system. 

The transformation of the velocity vector in the local coordinates system to the inertial 
frame may be accomplished using the following: 



cos (x*,x) cos (y«,x) cos (z*,x) 



= 

cos (x*,y) cos (ye,y) cos (z*,y) 


St 


cos (x*,h) cos (y*,h) cos (z*,h) . 


. 0. 


(A.9) 


The direction cosines can be obtained as 



r 


-* ■ 

' it ' 


A 

r-* 
>— * 

V 

A 

t3* 

V 

<n , i> 




— * — •* 

, 

y 

— 

< tl , j > 

A 

►s* 

t— * 

V 

< n , j > 

_ t . 


. <"ti , k> 

<"t 2 , k> 

<~n , k> . 


'yi 

o 


(A. 10) 


where < •, • > denotes the inner product. 

Thus, the kinematic equations of motion over the terrain profile becomes 

. Vcosx fJyVsinX 

Vl+f* Vl+fxVl + fJ+fy 


(A.11) 
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VTTff V sin x 

Vi+^+fJ 


(A.12) 


z = f(x,y) (A. 13) 

These equations were used in the trajectory planning schemes discussed in Chapters HI and 
IV. 
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APPENDIX B 

NUMERICAL CONJUGATE POINT TESTING 
FOR FIXED END-POINTS PROBLEMS 

B-1 Introduction 

The need for conjugate point test in the calculus of variations arise from the 
requirement that the second variation evaluated along an optimal path must be greater than 
zero for all admissible variations of states and costates. A pointwise test for the sign of 
second variation is the well known Legendre-Clebsch necessary condition. For extremals 
of finite length, however, the task of ensuring that the second variation is positive for all 
admissible neighboring paths is more involved [73]. The proposed Jacobi test [97] seeks 
the min imum value of the second variation. This problem is called the accessory minimum 
problem. The objective here is to find a system of variations which gives the value of zero 
to the second variation. If this is possible, it implies that a neighboring path is competitive 
and that the extremal furnishes at best an improper minimum of performance index and at 
worst merely a stationary value. 

Analytical conjugate point test is impossible in all but very simple optimal control 
problems. However, a family of numerical methods are available in the literature. This 
appendix examines three of these techniques. The first two are very general, while the 
third approach is useful in more classical situations such as the optimal trajectory planning 
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problem discussed in Chapter III. All these methods are equivalent, although some are 
numerically more efficient than others. 


B.2 Second Variation 

Consider the general optimal control problem described in Chapter II. 


* = f(x,U,t) , x(to) = Xo 

(B.l) 

i = -H x 

(B.2) 

II 

o 

(B.3) 

H (x,u,A.,t) = L (x,u,t) + X T f (x,u,t) 

(B.4) 

P(x(t f ),t f )=0 

(B.5) 

Q(x(t f ),t f ) = g(x(t f ),t f ) + v T P(x(t f ),t f ) 

(B.6) 


(B.7) 


The augmented performance index is given by: 

J a * Q(x(tf),tf) + f {H(x,u,A.,t) - X T x} dt 

J to 

Expanding the augmented performance index (B.8) to second order [12], one has 
AJ * = 5J * + ^ 5 2 J a + higher order terms 
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= Q5x(t f ) - X T (t f )5x(t f ) + f [(K+ H x ) 8x+H u 8u] dt + ^8x T( t f )|^8x(t f ) 




[V 


5u T 


H xx H xu 
H„ H^j J 


8x 

L5uJ 


dt + h.o.t. 


(B.9) 


Since the first variation terms vanish along every stationary trajectory, second-variation 
becomes 


8V = Sx^q^x^) + | [ 5x t 5u t] 

OX 1 


c 


^xx ^xu 


H U x H, 


uu 


5x 

-8u- 


dt 


(B.10) 


The differential constraints and boundary conditions in this problem are 


8x = f x 5x + fjiu. 

(B.ll) 

8x(0) = 0, 

(B.12) 

8P = [P x 8x] tf = 0 

(B.13) 


The accessory minimum problem attempts to find the minimum value of the second 
variation (B.10) subject to the differential constraints (B.ll) and boundary conditions 
(B.12), (B.13). The Euler-Lagrange equations for this accessary minimum problem are 
given by 

bi = - Hxx5x - f!5X. - H xu 5u (B . 1 4) 

H ux 5x + fj8X + H uu 5u = 0 (B. 1 5) 

8X(t f ) = {g„ + (v T P„)J8x| 

H 


(B.16) 
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If one assumes that control is nonsingular (i.e., H uu * 0), the expression (B.15) may be 
put in the form 


5u = -H u ^H ux 5x + fj5X) 


(B.17) 


In this case, one may substitute for 5u in (B. 14) and (B, 1 1) resulting in 


ffcl 


A(t) 

-B(t) 

5x 

1 

to 
1 


. -C(t) 

-A T (t) . 

-8X. 


(B.18) 


where. 


A(t) = f x -fJi^H ux , 


B(t) = fuHuifJ, 

This set of linear differential equations must satisfy the boundary conditions (B.12), 
(B.13), and (B.16). 


B.3 The Riccati Eauati oiLJVlethod 

* 

Since the accessary minimum problem (B.18) is linear and the differential equations 
and the terminal boundary conditions are homogeneous, both 5x(t) and 5X(t) are 
proportional to Sx(to), or proportional to 5X(tf). Using the backward sweep [41], where in 
one assumes a solution of the form 


8X(t) = S(t) 5x(t), S(tf) = Q* x 


(B.19) 
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a matrix Riccati equation can be obtained as follows [71]: 
§ = - SA -A t S + SBS - C, S(tf)=Qx X 


(B.20) 


Existence of a bounded symmetric-matrix S to the matrix Riccati equation (B.20) 
determines the existence of a conjugate point. In other words, if S becomes infinite at any 

A 

point t along a test extremal, then the second variation is zero in the interval [to. t]. In such 
a case, a neighboring extremal may furnish a lower value of the performance index. 

B.4 The Kellv-M over Method 

Regarding the linearized Euler equations (B. 18) as a mapping between the variations 
of unspecified initial costates, 8Xj(0) (i = 1,2, 3,—, n), and the variations of states, 8xj(t) 

(i = 1, 2, 3,—, n), the following can be written 


8xj(t) 



8^(0) 

8x2(0 


9 xt 8x) 9 xj 


8X 2 (0) 

: 


d^i o dX 2g dKo 

, • * • 

. • • ' » 


; 



dXn 





dx u ax 2o dKo _ 



_ 8x n (t) _ 



8X^(0) 


(B.21) 


where X.i 0 s Xi(0), i = 1 , 2 ,..., n 
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Setting 


9xi dxi 


axi n 


a^io 20 




M(t) = 


ax n ax n 

ax.io b ^ 2 o 


ax n 

a^nc _ 


(B.22) 


Kelly [73], Cicala [75], and Moyer [98] showed that the rank of the matrix M(t) 
provides the criterion for the existence of a conjugate point If the rank of the matrix M(t) 
drops at any point along the extremal, it indicates the existence of a conjugate point. In 
other words, neighboring extremals starting at to in n-dimensional state space must collapse 
into a smaller dimension at a conjugate point 


B.5 The Bliss Method 

In Reference 70, a numerical method for conjugate point testing was suggested. For a 
two dimensional system, if four fundamental solutions (x^yO (i = 1,2, 3,4) to the linearized 
Euler-Lagrange equations can be obtained, then the determinant of the matrix 


Xl 

X2 

X3 

X4 

yi 

Y2 

y3 

y4 

*i 

*2 

X3 

*4 

?i 

yi 

y3 

y4 


d(t) = 


(B.23) 
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will be non- zero. Since (xi,yi) (i = 1,2, 3, 4) are the solutions of the two dimensional 
accessory minimum problem, their linear combination is also a solution of the accessory 
system. Thus, 


x = ClXi + C2X2 + C3X3 + C4X4 

y = ciyi + C 2 y 2 + C 3 y3 + C4y4 


(B.24) 


This fact is exploited in the following theorem. In (B.23), note that (xj.yO (i - 1,2,3, 4) 
can be computed from the linearized costate equations. 


Theorem In a two dimensional variational problem, if four solutions (x L yi) (i — 
123 , 4 ) of the accessory equations of a nonsingular extremal arc are formed as a matrix 


and if the determinant of matrix 


D(t,t 0 ) = 


x^t) x 2 (t) 

y j(t) y 2 (t) 

Xi(t 0 ) x 2 (t 0 ) 

y i(t 0 > y 2 ( t o) 


x 3 (t) X 4 (t) ■ 

y 3 (t) y 4 (t) 

XjCto) x 4 (tg) 

y 3 (t 0 ) y 4 ( t o> 


(B.25) 


is not identically zero, then the conjugate point to point(x(to),y(to)) is the (x(t),y(t)) which 
makes determinant of matrix D(t, to) zero. 

Since the initial conditions on xi, X 2 , yi, y 2 arc arbitrary, matrix D(t,to)'can be re- 
expressed as 

|’x 1 (t) x 2 (t) x 3 (t) x 4 (t) ' 
y x (t) y 2 (t) y 3 (t) y 4 (t) 

0 0 x 3 (t Q ) x 4 (t Q ) 

0 0 y 3 (t 0 ) y 4 (t 0 ) 


D(t,t 0 ) = 


(B.26) 
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The determinant of this matrix D can be changed using Laplace Expansion as follows: 


A(t,to) = 


xi(t) x 2 (t) 
yi(t) y 2 (t) 


X3(to)X4(to) 

y3(to)y4(to) 


(B.27) 


The second determinant of right-hand side in above equation (B.27) is arbitrary. Hence, 
for A(t,to) to be zero, the characteristic determinant 

(B.28) 


A(t) = 


xi(t) 

yi(t) 


X2(0 

y2(t) 


should be zero at conjugate point t = t*, < t* < t^. 


M Relationship 

Jacobi's differential equations, or accessory system of differential equations, are 
linear and homogeneous. There exist two well-known methods for solving linear two- 
point boundary- value problems [12]: the backward sweep method and the transition 
matrix method. The backward sweep method results in the Riccati equations. On the 
other hand, Kelly-Moyer and Bliss Methods are based on the state transition matrix. 

Relationship of Riccati Equation Method and Keliv-Mover Method 

To test conjugate point by Kelly-Moyer method, the sign of determinant of matrix M(t) 
should be evaluated along the extremal. Equation (B.18) can be rewritten in terms of the 
partitioned transition matrix <I>(t,to) as 


8x(t) 


^xx(t»to) 

<*>xX(t,to) 


Sx(to) 

5X(t) . 


. G>Xx(t,to) 

^uCt.to) . 


. mo) . 


(B.29) 
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Using the expression (B.21), one can find out that 


M(t) = <D x x (Mo) 


(B.30) 


From the forward sweep equation, one has [12] 


5X(to) = S(to) 5x(to) 


(B.31) 


From this. 


M _1 (t) = S(t,to) 


(B.32) 


Relationship of Kellv-Moyer Method and Bliss Method 

The following corollary given by Bliss illustrates the relationship between the Kelly- 
Moyer method and Bliss method. 

Corollary In a two dimensional variational problem, if the accessory equations of 
a nonsingular extremal arc are a trajectory of a four-parameter family 

x(t,a,b,c,d), y(t,a,b,c,d) (B.33) 


then the conjugate point to a point to is determined by the zeros t of the determinant of 


matrix 


D(t,to) 


x a (t) x b (t) Xc(t) Xd(t) 

y a (t) yb(t) y c (0 yd(0 

X a (to) Xb(to) Xc(to) Xd(to) 

y a (to) yb(to) y c (to) yd(to) . 


(B.34) 


provided that the determinant of matrix D is not identically zero along trajectory. 


IL7 Numerical Effort 


All matrices mentioned in the foregoing are n by n. Thus, there is no difference in the 
size of matrix. To find the existence of conjugate point, however, the first method based on 
backward sweep checks the boundedness of a matrix, but the methods based on the 
transition matrix, Kelly-Moyer method and Bliss method, check the sign of determinant of 
the matrix. The latter approach is numerically better conditioned. 

The matrix Riccati method requires the integration of n x n equations and thus 
demands a formidable amount of calculations. Primary difficulty here is the requirement 
for various partial derivatives. In the present case, this would mean the computation of 
higher order terrain gradients. 

The numerical evaluation of the elements of Kelly-Moyer method evidently requires 

computer codes considerably more complex than those required for the calculation of Euler 

solutions. As an alternative to partial derivatives, Menon [99] used a scheme in which the 
partial derivatives with respect to the are calculated approximately in terms of difference 

quotients. Thus, small increments in the are employed in the evaluation of neighboring 
solutions of the extremal. 

The numerical Bliss method is similar to the Kelly-Moyer method. Instead of partial 
derivatives, fundamental solutions are obtained numerically. As an alternative to 
fundamental solutions, in the present work a computer code was assembled based on the 
finite differences. This code generates three trajectories corresponding to each heading 
angle, the first being the nominal and next two, the neighboring trajectories obtained by 
perturbing the initial value of the heading angle by a small amount in the positive and 
negative sense. This is equivalent to perturbing the initial values of costates while 
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enforcing the constant of motion. The required fundamental solutions are then computed 
using a finite difference scheme. 

The main difficulty encountered in both finite difference quotients computations [99] 
and finite difference method is the errors arising from higher order effects. These errors 
can be controlled to a certain extent by verifying the linearity of the x(t) and y(t) differences 
versus the magnitude of the corresponding increment in the initial heading angle. This 

check can be incorporated in the computer program. 

If the final end-point is not fixed, the matrix Riccati method requires additional matrix 

calculation, see Reference 12 for further details. 
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APPENDIX C 


SEPARABILITY OF THE HAMILTONIAN AND ITS 
CONSEQUENCE ON DIFFERENTIAL GAME SOLUTIONS 


C-1 Introduction 

While discussing the trajectory planning schemes in Chapter IV, it was pointed out that 
in order to ensure that the outcome of the game is not influenced by the order of action, the 
value of the game should exist. The value exists only when the final performance index 
has a saddle point. The conditions for the existence of the saddle point are the separability 
of the variational Hamiltonian and continuous mapping. In the following, each of these 
issues will be elaborated. 


C.2 Separability of The Hamiltonian 

Consider a system of differential equations written in vector form 

x = f (x,<(>,Y,t), x(t 0 )=x 0 

where, x(t) : = state vector of dimension n, x e X 

<|>(t) : = control of Player 1 of dimension 4, <t> e <X> 
y(t) : = control of Player 2 of dimension m, \|f e 'P 


with terminal constraints: 



P (x(tf),t f ) = 0, t f is free 


(C.2) 


and the performance index being defined as 

( U 

J[x,<J>,V|/,t] = g(x(tf),tf ) + I L(x,<t>,\|/,t) dt (C.3) 

Jio 

The variational Hamiltonian for this problem is given by 

H = L + X T f + pi T P + \i2 g (C.4) 

If the right-hand side of system dynamics (C.1) and the integral part of payoff (C.3) are in 


the form : * "" r " .. •“ 

f (x,4»,v,t) = fj (x,4>,t) + f 2 (x,\jr,t) (C.5) 

L (x,<t»,v,t) = Lj (x,<J>,t) + (x,\y, t) (C.6) 

then the Hamiltonian is separable in the space of controls \|r and <|> as follows: 

H = H 1 + H 2 + p 1 T P + p 2 T g (C.7) 

Hj = Lj (x,<|),t) + X T fj (x,4>,t) , (C.8) 

H 2 = L 2 (x,y,t) + X T f 2 (x,\j/,t) (C.9) 


In this case, Friedman [39] showed that there exists a saddle point for the Hamiltonian. 
The order of action is immaterial in this case, i.e., 

max min H(x,<f>,\jr, A,t) = min max H(x,<|>,V,X,t) = H(x,6°,\ir°,X,t) 

V 4> 4> V 

(C.10) 
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Equation (C.10) is known as the Isaacs condition in differential game [39]. This condition 
is also known as the saddle point condition. 


C .3 Sufficient Condition fo r a Game to have Value 

As discussed in the previous section, the Isaacs condition guarantees the existence of 
the pointwise saddle point. If the final performance index in the game has a saddle point, 
then one is guaranteed that the order of action is immaterial. This section treats how the 
pointwise saddle point can be transformed to the global saddle point. The chief 
requirement here is that the game must have a value. 

Friedman [39] [100] proved that a sufficient condition for the game to have a value is 

the Isaacs condition and the following requirements.: 

(i) f(x,<f>,y,t) is continuous and satisfies the Lipschitz condition, 

(ii) g(x(t f ), tj) and L(x,<t>,y,t) are continuous. 

Here, condition (i) is required so that the trajectory generated by integrating equations of 
motion (C.l) is unique and continuous, while condition (ii) is required to ensure that the 
performance index is continuous. Under these conditions, if the variational Hamiltonian is 
separable, the performance index will be a unique quantity regardless of the player's order 
of action. 

To show the existence of the value, Friedman [39] and Fleming [101] considered an 
approximating upper and lower game for every partition of the time interval [to, tf] into K 
equal sub-intervals. Using the requirement (i), the trajectory evolves according to 
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(C.11) 

Two possibilities are next considered. The first one in which the minimizing player 
plays first. In the second case, the maximizing player plays first. The performance index 
in each case is 

V K + = inf sup • • • inf sup J [x, <(> ., y.,t] j =0, 1,2,- • ,K- 1 (C.12) 

w. A 6 «, J J 

Y# 4>, 

V K ~ = sup inf • - sup inf J[x,4>,\jr .t] j =0,1,2,- • -,K- 1 (C.13) 


J ftj+i 

f(x,<(>j,\Pj,T)dT 
< j 


With the increment of K, 

lim V+ = V + , lim V2 = V 
K-*» K-»~ K 


It can be shown [39] that in the general case, 
v + 2> y- 


(C.14) 


(CAS) 


Note that the above expression will be an equality if the Isaac's condition is satisfied. In 
such a case. 


APPENDIX D 


CONVEXITY CONDITION FOR GLOBAL MINIMUM 


Convex functionals play a special role in optimal control theory because most of the 
theory of local extrema for general nonlinear functionals can be strengthened to become 
global theory when applied to convex functionals [69]. 

By definition in Reference 35, a real-valued functional J: Y — » 9f is convex if 

J[a yi + (1 - tx) y 2 ] £ a J[yi] + (1 - a) J[y 2 ]» V yi,y 2 e 7 , 0 < a < 1 (D.l) 

The following theorem [35] asserts the global nature of results for minimization problems 
with convex functionals such as those defined in (D.l). 

Theorem If J(yo) is a weak local minimum of J(y) on Y and J is convex at yo relative 
to Y, then J(yo) is a global minimum. 

A sufficient condition for convexity of functionals is that the integrand be convex for 
the integral interval. Thus the test of convexity for a functional can be reduced to a test of 
the convexity of the integrand. Typical integrands that are convex appear in Linear 
Quadratic Regulator theory and in the accessory minimum problem related to the Jacobi 
condition in nonsingular cases. 
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Table 3.1 Digital Terrain Data used for Trajectory Planning 


E 


121 


ON 


00 


NO 


m 


co 


oo 

CS 

Tf 


8 

CO 


CO 


© o o 

OO M (N 
CS 


8888 

(N (N CO 


r- 

5 


— • o o 

M N 
VO h* O 


o 

m 

r- 


cs 

Tf _ 
CS Tf 


8 


ON 

00 

o 


,-4 i-4 *— I i— < i— « i-H i-H i-H i-H 1— < i“H i-H i-H w—t CS 


8 8 S 8 

- -- — 


CO CO 


o o 

CS CS 


o 

Tf 

CS 


O 

tj- 

cs 


5 

cs 


8 


o 

00 

cs 


CS 


CS 


© © © 

cs cs o 

H PK (S 


a 


CS *-» 


o q 
cs 


cs 


P © 


s 


o 

cs 


I 


_ r- 
3 O vo 
O CS CS 
• • • 
^ 

883 

CS CS ^ 


o o o 
cs cs cs 


8 0 CO CO 

oo in no 

cs CS CO CO 


o o © cs 
TT Tf Tf NO 
CS CS <S CO 


O CS o ^ 
CS 1— I Tt VO 
^ CS (S CO 


OOOO 

CS Tf Tf 25 

1-1 cs CS CO 


o O © © 

cs Tf Tf 25 

1-t CS CS CO 


OOOO 

saa? 


cs 

NO 

CO 


3. 

I 


cs 

cs 

m 


cs 

5 


8 

m 


o 

cs 

m 


o o 
00 00 


CO CS O 
NO OO ^ 

r- no no 


88 

r- cs 

^ cs 

CS O 

no m 


cs 

8 


CS 

00 CS Q 

no m vo 


Tf 

00 

Ti- 


es cs cs 

^ H 

O CS NO 

Tf CO 1— J 

cs cs cs 

1— ' CO O 
CS Tf cs 
CO CS CO 

cs cs cs cs’ 

CO rH Q 
CS Ti ® 
CO CS i— 1 


CS on o 

NO OO CS 

m tj- m 


On m O 
oo oo oo 

Tf Tf Tf 


m cs co 
00 oo 00 

Tf Tf Tf 


B 8 

r- *-» 

1 -H cs 

CO O 
VO OO 

r- *-> 

• • 

*-< cs 

CS ON 

NO © 

q 

i-4 cs 

^H 

00 Tf 

VO 00 


On 


CS CS CS 

cs ^ © 
Tt 00 cs 
cs o * 
cs cs 
cs cs 

Tf NO 

oo r- 


cs cs 
cs cs 
m m 


ooooooooo 

CSl TfTfCSCSCSTfOOOOCS,,- 

OOHHH(S(Scj»n»n^vo 3 


CQNCS^CS^^ON-; 


m Tf 


m »n 


<u 

.D 


<ul 

01 

c 

03 


O On On 
00 CS CS 
On h - 1 

o ^ ^ 

o q q 
CS cs cs 


88888 

cs cs cs co m 


q o 
cs ^ 


m 


co 

8 


OOOOOOO^ 
TfTfTfTtOOOOCSCS 
CS CS CS ^T 


r- *n 
O o 

Tt Tf 


8 

no no m m Tf 


s 

Tf 


i— • Tf 

CS cs 
in m 


CO OO i-H 

cs cs cs 

in m in 


cs «-r o 

3 $> § 

r— < 

i— < i— < CS 

8 2 3 


c 

o 


c§ 


o 

On | 2 


OO 


8 

cs 


88 ^ 

CS CS H 


o © 

m Tf 
<S cs 


88 

cs co 


S O oo 
© cs 
co 05 oo 


o o 

cs OO 
wo VO 


g o <0 (N H H 


_ CS O O O CS 

t^mTfTfTfminwo 


o q 
cs 


r- O o ON cs 
CS CS ^ oo oo 


cs cs in 


Tf cs O © oo Q CS 
Tf uo Tf Tf Tf no m 


NO NO 


CO O CS CS no 

cs -h OO 00 NO 

in no no no 


NO 


3 


© 

CS 


© 

m 


cs 


© 

00 

cs 


It 
i s 

co co 


© © © 

cs cs © 

m m Tf 


88 ? 

co co co 


i— i © © © co On i-« 

© © OO cs cs ^ cs 

Tf Tf ^ m no no *n 

» « i * » • • 

i-H i-H i— I 1— » i — 1 i— ( i-H 

CS © © ON CS *-« © 

© 00 CS oo CS NO Q 

Tf Tf m On On r- vO 


3 3 S ? S 

25 cFn oo r^ 


Tf 1— © CS o 

cs OO 00 00 00 

VO vo vo vo NO 


m 


cs 


© 

m 

cs 


as 

Tf CO 


8 ©©m©©o^i— 'O 
Tfrf©oooooooo ^3 
Tfcncn^Tfvoqqqq 


vo © CN 

Tf cs © 00 OO 

vo r- r- no no 


© 

m 

cs 


© 

00 

Tf 


© © 
CS Tf 

m co 


© © © 

Tf Tf Tf 

co co co 


80 S 

Tf Tf Tf 


1—1 1— 1 o 
00 CS NO 

Tf no r- 


1 — < CS 1 — ' © cs 

NO VO Tf Tf Tf 

r- 00 00 00 00 


co 


cs 


co co 


cs co 


© 

00 

cs 


s 

CO 


CO 


I 

« 

s 

CO 


© 

00 

cs 


9 

CO 


© cs 
CS co 
m r- 


888 

CO CO Tf 


m © © © 

© 00 00 © 

Tf Tf Tf ® 


§ © 

. ^ 

r-H i-H 

ss 

CO CO 

ii 


888 

CO CO Tf 


in 


So s a 

NO r-' On 


8 

r- 


© *-H © cs cs 

^ cs ^ cs cs 

00 On © © ON 

^4 i -4 cs’ cs i -4 


888 
Tf Tf Tf 


88 88 

in ® h o 

1 -H rH *— « CS 

8S88 

vo © m 
• p • • 
1-H 1-H CS CS 


© i-H On 

£ ^0 © 
© ^ © 
CS CS CS 

CO 1 -H 1— ■ 

S cs 00 

Tf © 

cs cs cs 

cs q 
cs cs cs 


^ Ov © 

8 Tf v5 

CO ^ 

cs cs cs 
cs cs © 


CO 

cs 

ON 


cs m 


vb cs 60 Tf 

i-H On OO 00 

CS CS 1-4 1-4 r-4 


cs 

Tf 

r- 


5 ? 


CS Tf 


00 00 00 00 



aSUE^-UMOQ 




122 


Table 5.1 Maximum Rate of Climb (ft/sec) for AH- IS 


Speed (ft/sec) Density Altitude 



sea level 
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2000 ft 
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Figure 3. 1 The Coordinate System 
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Figure 32 Flow Chart for Generating Euler Solutions 












Figure 3.3 Sample Terrain Map of the Nassau -Valley, California 
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Figure 3.4 Euler Solutions for Minimum Flight Time Criterion (K = 0) 
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Figure 3.6 Comparison between Minimum Time Trajectory 
and Maximum Terrain Masking Trajectory 
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Figure 3.7 Altitude Profiles along Different Criterion Trajectories 
Solid lme:maximum terrain masking trajectory 
Dotted line:minimum flight time trajectory 
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Figure 3.8 Characteristic Determinant A(t) along the Trajectory A 
for Maximum Terrain Masking (K = 0.99) 
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Figure 3. 10 Comparison of Performance Index along Trajectories A and B 
for Maximum Terrain Masking 



Figure 3. 1 1 Two Extremals for Minimum Flight Time Criterion 
having a Same Terminal Position 
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Figure 3.14 Euler Solutions for Optimal Route Planning No.2 
(<x= 10**5, e =0.001) 




Figure 3.15 Comparison between Straight Trajectory and an Optimal 
Route Planning No.2 Trajectory (a = 10**6, e =0.01) 
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Figure 3.16 Altitude Profiles along Straight Trajectory and Optimal Trajectory 
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Wind 



Considering Wind Effects (u/V =0.1) 
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Wind 



Figure 3.18 Euler Solutions for Maximum Terrain Masking Criterion 
Considering Wind Effects (u/V =0.1) 
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Figure 4. 1 Trajectories for the Pursuer and Evader 
(W P = 0.0,Wc = 0.0) 
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Figure 4.2 Trajectories for the Pursuer and Evader 
(W P = 1.0, We = 1.0) 
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Figure 4.3 Altitude Histories for the Pursuer and Evader 
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Figure 4.4 Trajectories for the Pursuer and Evader 
(Wp = 1 .0, We = 1 .0) 
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Figure 4.6 Trajectories for the Pursuer and Evader 
(W P = 0.0, W c = 1.0) 
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Figure 4.7 Trajectories for the Pursuer and Evader 
(W P = 0.5, We = 1.0) 
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Figure 4.8 Legendre-Clebsch Test for the Pursuer and Evader 
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Figur 4.9 Legendre-Clebsch Test for the Pursuer and Evader 







Figure 4. 1 1 Trajectories for the Pursuer and Evader 
(c = 1.0, a = 0.0007, (5 = 0.00001) 
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Figure 4.12 Speed Histories 


for the Pursuer and Evader 
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Figure 4. 13 Load Factor Histories for the Pursuer and Evader 
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Figure 4.14 Bank Angle Histories for the Pursuer and Evader 
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Figure 4.17 Trajectories for the Pursuer and Evader 
(o = 1.0, a = 0.0007, P = 0.00001) 
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Figure 4.18 Speed Histories for the Pursuer and Evader 
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Figure 4.19 Load Factor Histories for the Pursuer and Evader 
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Figure 4.20 Roll Attitude Histories for the Pursuer and Evader 
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Figure 4.21 Pitch Attitude Histories fot the Pursuer and Evader 
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Figure 5.2 Longitudinal Cyclic Control for Maximum Masking Trajectory 
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Figure 5.3 Lateral Cycli< 
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Figure 5.6 Pitch Attitude Response for Maximum Masking Trajectory 
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Figure 5. 10 Longitudinal Cyclic Control for Minimum Time Trajectory 
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Figure 5. 1 1 Lateral Cyclic Control for Minimum Time Trajectory 
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Figure 5. 13 Collective Control for Minimum Time Trajectory 
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Figure 5.14 Pitch Attitude Response for Minimum Time Trajectory 
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Figure 5.15 Roll Attitude Response for Minimum Time Trajectory 



Rate of Climb (ft/sec 


r 


171 



Time (sec) 

Figure 5.9 Rate of Climb for Minimum Time Trajectory 
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